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