1 subroutine grad1_so 2 $ ( H, lbuf, scr, lscr, dens, wdens, frc_nuc, 3 $ frc_kin, frc_wgh, 4 $ g_dens, g_wdens, basis, geom, nproc, nat, 5 $ max_at_bf, oskel, 6 & frc_so, densx, densy, densz) 7c$Id$ 8 9C one electron contribution to RHF, ROHF and UHF gradients 10C now also UMP2 11 12 implicit none 13 14#include "mafdecls.fh" 15#include "global.fh" 16#include "geom.fh" 17#include "bas.fh" 18#include "rtdb.fh" 19#include "sym.fh" 20 21C-------------------------parameters-------------------------------- 22 integer lbuf, lscr, 23 $ g_dens(6), ! density matrix (summed if ROHF, UHF) 24 $ g_wdens, ! weighted density (Lagrangian) 25 $ basis, geom, nproc, nat, max_at_bf 26 27 double precision H, ! integral derivatives 28 $ scr, 29 $ dens, densx, densy, densz, ! local density block 30 $ wdens, ! local weighted density block 31 $ frc_nuc, frc_kin, frc_wgh, frc_so ! forces arrays 32 33 dimension H ( lbuf ), frc_nuc(3, nat), frc_kin(3, nat), 34 $ frc_wgh(3, nat), frc_so(3, nat), scr(lscr), 35 $ dens(max_at_bf,max_at_bf), wdens(max_at_bf,max_at_bf), 36 $ densx(max_at_bf,max_at_bf), densy(max_at_bf,max_at_bf), 37 $ densz(max_at_bf,max_at_bf) 38 39 logical oskel ! symmetry? 40 41C-------------------------local variables-------------------------- 42 43 integer ijatom, next, iat1, iat2, iat3, ish1, ish2, 44 $ iab1f, iab1l, iab2f, iab2l, iac1f, iac1l, iac2f, iac2l, 45 $ if1, il1, if2, il2, 46 $ icart, ic, nint, ip1, ip2 47 48 double precision crd1, crd2 ! atomic coordinates 49 dimension crd1(3), crd2(3) 50 51 integer idatom 52 dimension idatom(2) 53 54 double precision dE, dx, dy, dz, qfac, fact, q1, q2 55 56 logical status, pointforce 57 58 character*16 name 59 integer nxtask, task_size 60 external nxtask 61 62 task_size = 1 63 status = rtdb_parallel(.true.) ! Broadcast reads to all processes 64 65 pointforce = geom_include_bqbq(geom) 66 67 call hf_print_set(1) 68 69 ijatom = -1 70 next = nxtask(nproc,task_size) 71 do 90, iat1 = 1, nat 72 do 80, iat2 = 1, iat1 73 74 ijatom = ijatom + 1 75 if ( ijatom .eq. next ) then 76 77 status = bas_ce2bfr(basis,iat1,iab1f,iab1l) 78 status = bas_ce2bfr(basis,iat2,iab2f,iab2l) 79 80 if (iab1f.le.0 .or. iab2f.le.0) then 81c 82c At least one center has no functions on it ... next atom 83c 84 goto 1010 85 endif 86 87 if (oskel) then 88 if (.not. sym_atom_pair(geom, iat1, iat2, qfac)) 89 $ goto 1010 90 else 91 qfac = 1.0d0 92 endif 93 94 status = bas_ce2cnr(basis,iat1,iac1f,iac1l) 95 status = bas_ce2cnr(basis,iat2,iac2f,iac2l) 96 97 call ga_get (g_dens(1), 98 & iab1f,iab1l,iab2f,iab2l,dens,max_at_bf) 99 call ga_get(g_wdens, 100 & iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf) 101 call ga_get (g_dens(3), 102 & iab1f,iab1l,iab2f,iab2l,densz,max_at_bf) 103 call ga_get (g_dens(4), 104 & iab1f,iab1l,iab2f,iab2l,densy,max_at_bf) 105 call ga_get (g_dens(5), 106 & iab1f,iab1l,iab2f,iab2l,densx,max_at_bf) 107 do 70, ish1 = iac1f, iac1l 108 if ( iat1.eq.iat2 ) iac2l = ish1 109 do 60, ish2 = iac2f, iac2l 110 111C shell block in atomic (D/Dw)-matrix block 112 status = bas_cn2bfr(basis,ish1,if1,il1) 113 if1 = if1 - iab1f + 1 114 il1 = il1 - iab1f + 1 115 status = bas_cn2bfr(basis,ish2,if2,il2) 116 if2 = if2 - iab2f + 1 117 il2 = il2 - iab2f + 1 118 119 nint = ( il1 - if1 + 1 ) * ( il2 - if2 + 1 ) 120 121C overlap derivatives 122 call intd_1eov(basis,ish1,basis,ish2,lscr,scr, 123 & lbuf,H,idatom) 124 125C Dw x S 126 127 if ( idatom(1) .ge. 1 ) then 128C idatom(1).ge.0 <=> idatom(2).ge.0 (no check necessary) 129 ic = 1 130 do 28, icart = 1, 3 131 de = 0.D0 132 do 22, ip1 = if1, il1 133 do 20, ip2 = if2, il2 134 dE = dE + wdens(ip1,ip2) * H(ic) 135 ic = ic + 1 136 20 continue 137 22 continue 138 dE = dE * qfac 139 frc_wgh(icart,idatom(1)) = frc_wgh(icart,idatom(1)) 140 $ - dE - dE 141 frc_wgh(icart,idatom(2)) = frc_wgh(icart,idatom(2)) 142 $ + dE + dE 143 28 continue 144 endif 145 146C 1el. derivatives 147 call intd_1eh1(basis,ish1,basis,ish2,lscr,scr, 148 & lbuf,H) 149 150C D x H 151 152 ic=1 153 do 50, iat3 = 1, nat 154 do 40, icart = 1, 3 155 dE = 0.D0 156 do 31, ip1 = if1, il1 157 do 30, ip2 = if2, il2 158 dE = dE + dens(ip1,ip2) * H(ic) 159 ic = ic + 1 160 30 continue 161 31 continue 162 if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE 163 dE = dE * qfac 164 frc_kin(icart,iat3) = frc_kin(icart,iat3) + dE 165 40 continue 166 50 continue 167 168C 1el. so. derivatives 169 call intd_1eso(basis,ish1,basis,ish2,lscr,scr, 170 & lbuf,H) 171C Dso x Hso 172 173 ic=1 174 do 150, iat3 = 1, nat 175 do 140, icart = 1, 3 176c z componet 177 dE = 0.D0 178 do 131, ip1 = if1, il1 179 do 128, ip2 = if2, il2 180 dE = dE - densz(ip1,ip2)*H(ic)*0.5d0 181 ic = ic + 1 182 128 continue 183 131 continue 184 if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE 185 dE = dE * qfac 186 frc_so(icart,iat3) = frc_so(icart,iat3) + dE 187c y componet 188 dE = 0.D0 189 do 230, ip1 = if1, il1 190 do 231, ip2 = if2, il2 191 dE = dE - densy(ip1,ip2)*H(ic)*0.5d0 192 ic = ic + 1 193 231 continue 194 230 continue 195 if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE 196 dE = dE * qfac 197 frc_so(icart,iat3) = frc_so(icart,iat3) + dE 198c x component 199 dE = 0.D0 200 do 250, ip1 = if1, il1 201 do 251, ip2 = if2, il2 202 dE = dE - densx(ip1,ip2)*H(ic)*0.5d0 203 ic = ic + 1 204 251 continue 205 250 continue 206 if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE 207 dE = dE * qfac 208 frc_so(icart,iat3) = frc_so(icart,iat3) + dE 209 140 continue 210 150 continue 211 212 60 continue 213 70 continue 214 215 1010 continue 216 217C Vnn 218 219 if ( iat1 .NE. iat2 ) then 220 if (iab1f.ne.0 .or. iab2f.ne.0 .or. pointforce ) then 221C no forces between point charges (for John Nicholas) 222 status = geom_cent_get (geom, iat1, name, crd1, q1) 223 status = geom_cent_get (geom, iat2, name, crd2, q2) 224 dx = crd2(1) - crd1(1) 225 dy = crd2(2) - crd1(2) 226 dz = crd2(3) - crd1(3) 227 fact = q1 * q2 / SQRT ( dx*dx + dy*dy + dz*dz ) **3 228 dE = dx * fact 229 frc_nuc(1,iat1) = frc_nuc(1,iat1) + dE 230 frc_nuc(1,iat2) = frc_nuc(1,iat2) - dE 231 dE = dy * fact 232 frc_nuc(2,iat1) = frc_nuc(2,iat1) + dE 233 frc_nuc(2,iat2) = frc_nuc(2,iat2) - dE 234 dE = dz * fact 235 frc_nuc(3,iat1) = frc_nuc(3,iat1) + dE 236 frc_nuc(3,iat2) = frc_nuc(3,iat2) - dE 237 endif 238 endif 239 240 next = nxtask(nproc,task_size) 241 endif 242 243 80 continue 244 90 continue 245 next = nxtask(-nproc,task_size) 246c write(*,'("forces",9f10.5)')((frc_nuc(i,j),i=1,3),j=1,nat) 247c write(*,'("forces",9f10.5)')((frc_kin(i,j),i=1,3),j=1,nat) 248c write(*,'("forces",9f10.5)')((frc_wgh(i,j),i=1,3),j=1,nat) 249c write(*,'("forces",9f10.5)')((frc_so(i,j),i=1,3),j=1,nat) 250 251 return 252 end 253