1
subroutine get_vecF1(g_vecF1, ! out:
5
& vectors, ! in : MO vectors
6
& nbf, ! in : nr. basis functions
9
& npol, ! in : nr. polarizations
10
& lifetime, ! in : = (.true.,.false.) with/out damping
11
& nocc, ! in : nr. occupied MOs
12
& nvirt, ! in : nr. virtual MOs
13
& debug) ! in : = .true. for debugging
15
c Author : Fredy W. Aquino
17
c Note.- Modified from original aoresponse source code
18
c for extension to spin-unrestricted case
19
c original aoresponse source code was written by
20
c J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
22
c --> Experimental (not published yet)
27
#include "mafdecls.fh"
30
integer ispin,npol,ncomp,nocc(npol),nvirt(npol)
31
c Note.- g_vecF1(npol,ncomp) does not work for two dimensions
32
integer g_vecF1(2,2),g_vecF1_im(2,2),
33
& g_rhs(ncomp),g_rhs_im(ncomp),
36
integer ipm,ifld,ndir,nbf,nmo,disp
37
logical lifetime,debug
38
integer alo(3), ahi(3),
44
ndir=3 ! = nr directions (x,y,z)
45
c ---- Create (g_vecF1,g_vecF1_im) -------- START
54
write (cstemp,'(a,i1)') 'aor vecF1 ',ipm
55
if (.not.nga_create(MT_DBL,3,ahi,trim(cstemp),
56
& alo,g_vecF1(ispin,ipm) ))
57
& call errquit('fiao_f1: nga_create failed vecF1',0,GA_ERR)
58
call ga_zero(g_vecF1(ispin,ipm))
60
write (cstemp,'(a,i1)') 'aor vecF1_Im ',ipm
61
if (.not.nga_create(MT_DBL,3,ahi,trim(cstemp),
62
& alo,g_vecF1_im(ispin,ipm) ))
63
& call errquit('fiao_f1: nga_create failed E1_im',0,GA_ERR)
64
call ga_zero(g_vecF1_im(ispin,ipm))
67
enddo ! end-loop-ispin
68
c ---- Create (g_vecF1,g_vecF1_im) -------- END
77
if (.not.nga_create(MT_DBL,3,ahi,'U matrix',alo,g_u)) call
78
& errquit('giao_b1: nga_create failed g_u',0,GA_ERR)
96
disp=nocc(1)*nvirt(1)*(ispin-1)
98
phi(1) = disp+nocc(ispin)*nvirt(ispin)
99
qlo(1) = nocc(ispin)+1
110
c ======== Including g_u_ov (g_rhs --> g_u) ==== START
114
call nga_copy_patch('n',g_rhs(ipm),plo,phi,
116
c ======== Including g_u_ov (g_rhs --> g_u) ==== END
117
call nga_matmul_patch('n','n',1.0d0,0.0d0,
118
& vectors(ispin) ,blo,bhi,
120
& g_vecF1(ispin,ipm),clo,chi)
122
c ======== Including g_u_ov_im (g_rhs_im --> g_u_im) == START
124
call nga_copy_patch('n',g_rhs_im(ipm),plo,phi,
126
c ======== Including g_u_ov_im (g_rhs_im --> g_u_im) == END
127
call nga_matmul_patch('n','n',1.0d0,0.0d0,
128
& vectors(ispin) ,blo,bhi,
130
& g_vecF1_im(ispin,ipm),clo,chi)
131
end if ! end-if-lifetime (damping)
133
enddo ! end-loop-ifld
134
if (.not.ga_destroy(g_u)) call
135
& errquit('fiao_b1: ga_destroy failed g_d1',0,GA_ERR)
136
enddo ! end-loop-ispin
140
subroutine update_rhs_dipole(
142
& vectors, ! in : MO vectors
143
& rtdb, ! in : rtdb handle
144
& basis, ! in : basis handle
145
& lvelocity,! in : logical var
146
& nat, ! in : nr. atoms
147
& npol, ! in : nr of polarizations
148
& nocc, ! in : nr. occ shells
149
& nvirt, ! in : nr. virt shells
150
& nbf, ! in : nr. basis functions
151
& nmo, ! in : nr. MOs
152
& ncomp, ! in : nr components of ...
153
& debug) ! in : logical for debugging
155
c Author : Fredy W. Aquino
157
c Note.- Modified from original aoresponse source code
158
c for extension to spin-unrestricted case
159
c original aoresponse source code was written by
160
c J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
162
c --> Experimental (not published yet)
165
#include "errquit.fh"
167
#include "mafdecls.fh"
175
integer g_dipole,g_dipole1,g_rhs(ncomp)
177
integer vectors(npol)
178
integer ifld,ndir,nbf,nmo,ipm,nat,
179
& nocc(npol),nvirt(npol)
180
double precision origin(3)
181
data origin/0d0,0d0,0d0/
182
integer shift,disp,ispin
183
integer alo(3), ahi(3),
185
logical lvelocity,oskel,debug
188
ndir=3 ! nr directions (x,y,z)
189
c Get dipole integrals in GA and transform to MO set (virt,occ)
196
if (.not.nga_create(MT_DBL,3,ahi,'dip matrix',alo,g_dipole1)) call
197
& errquit('fiao_f1: nga_create failed g_dipole',0,GA_ERR)
199
if (.not.nga_create(MT_DBL,3,ahi,'dip matrix',alo,g_dipole)) call
200
& errquit('fiao_f1: nga_create failed g_dipole',0,GA_ERR)
201
if (debug) write (luout,*) 'fiao_f1: dipole matrix allocated'
202
c Get H10 in GA, using g_dipole array
203
c note: origin has been set to (0,0,0) for multipole integs.
204
call ga_zero(g_dipole1)
206
call int_giao_1ega(basis,basis,
207
& g_dipole1,'velocity', ! out
209
call ga_scale (g_dipole1, -1d0)
211
call int_mpole_1ega(basis,basis,
212
& g_dipole1,'dipole', ! out
215
if (debug) write (luout,*) 'fiao_f1: AO integrals done'
216
c ga_rhs(a,i) = ga_rhs(a,i) + H10(a,i)
217
c Transform H10 to MO and add to g_rhs.
218
c (the rhs is NOT divided by (e_a - e_i -/+ omega), this
219
c will be considered in the CPKS solver, in the precon-
220
c ditioner and the 1e part of the "product" routine)
221
c --------- g_dipole1 --> g_dipole --------- START
222
call ga_zero(g_dipole)
237
call nga_copy_patch('n',g_dipole1,blo,bhi,
239
enddo ! end-loop-ispin
240
c --------- g_dipole1 --> g_dipole --------- END
241
call giao_aotomo(g_dipole,vectors,nocc,nvirt,npol,ndir,nbf)
244
alo(1) = nocc(ispin)+1
250
disp=nocc(1)*nvirt(1)*(ispin-1)
252
bhi(1) = disp+nocc(ispin)*nvirt(ispin)
258
if (ga_nodeid().eq.0)
259
& write(*,*) '---- g_dipole-nw-------- START'
260
call ga_print(g_dipole)
261
if (ga_nodeid().eq.0)
262
& write(*,*) '---- g_dipole-nw-------- END'
265
call nga_add_patch(1.0d0,g_rhs(ipm),blo,bhi,
266
& 1.0d0,g_dipole ,alo,ahi,
267
& g_rhs(ipm),blo,bhi)
269
enddo ! end-loop-ispin
270
if (debug) write (luout,*) 'fiao_f1: dipole added to rhs'
271
c Cleanup g_dipole as we do not need it right now
272
if (.not.ga_destroy(g_dipole)) call
273
& errquit('fiao_f1: ga_destroy failed g_dipole',0,GA_ERR)
274
if (.not.ga_destroy(g_dipole1)) call
275
& errquit('fiao_f1: ga_destroy failed g_dipole',0,GA_ERR)
279
subroutine get_nocc(rtdb, ! in : rtdb handle
280
& nocc, ! out: nr occupations
281
& npol, ! out: nr of polarization
282
& nclosed,! in : nr closed shells
283
& nopen, ! in : nr open shells
284
& nvirt, ! in : nr virtual MOs
285
& scftyp, ! in : string = UHF or RHF
286
& ntot) ! out: sum_{i,npol} nocc(i)*nvirt(i)
288
c Author : Fredy W. Aquino
290
c Note.- Modified from original aoresponse source code
291
c for extension to spin-unrestricted case
292
c original aoresponse source code was written by
293
c J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
295
c --> Experimental (not published yet)
298
#include "errquit.fh"
299
#include "mafdecls.fh"
306
& npol,nocc(2),nclosed(2),
307
& nopen(2),nvirt(2),ntot
308
if (scftyp .eq. 'UHF') then
312
c ------ Store nopen in rtdb so that CPHF routine is happy ---- START
313
c In scf_init(): to avoid error message:
314
c ===> scf: no. of closed-shell electrons is not even!
315
if (.not. rtdb_put(rtdb, 'scf:nopen',
316
& MT_INT, 1, nocc(1)-nocc(2)))
317
* call errquit('get_nocc:rtdbput nopen failed',
320
c ------ Store nopen in rtdb so that CPHF routine is happy ---- END
321
else if (scftyp .eq. 'RHF') then
328
ntot=ntot+nocc(ispin)*nvirt(ispin)