1      subroutine hnd_sdfc(rtdb,geom,basis,vectors,nclosed,nopen,nvirt,
2     &                    nbf,nmo,pairlist,translate,ipairs,unique,
3     &                    i_pert,i_resp,tensor,coords,nfc,nsd,nsdfc)
4c $Id$
5      implicit none
6#include "errquit.fh"
7#include "global.fh"
8#include "mafdecls.fh"
9#include "msgids.fh"
10#include "geom.fh"
11#include "rtdb.fh"
12#include "bas.fh"
13#include "nwc_const.fh"
14#include "stdio.fh"
15#include "apiP.fh"
16c
17      integer rtdb                             ! [input] rtdb handle
18      integer basis                            ! [input] basis handle
19      integer geom                             ! [input] geometry handle
20      integer vectors(2)                       ! [input] vectors
21      integer nclosed(2), nopen(2), nvirt(2)   ! [input] occupation info
22      integer nbf, nmo                         ! [input] orbital info
23      integer ipairs                           ! [input] number of spin-spin pairs
24      integer pairlist(2*ipairs)               ! [input] list of the pairs
25      integer translate(2*ipairs)              ! [input] translation to unique list
26      integer i_pert, i_resp                   ! [input] # of unique responding and perturbing atoms
27      integer unique(i_pert+i_resp)            ! [input] list of unique atoms
28      double precision tensor(3,3,5,ipairs)    ! [output] spin-spin tensor, one for each term
29      double precision coords(3,i_pert+i_resp) ! [input] coordinates of unique atoms
30      double precision nfc, nsd, nsdfc         ! [input] prefactors for the three terms
31c
32      integer ixy, ix, iy, iz, ifld
33      integer alo(3), ahi(3), blo(3), bhi(3), clo(3), chi(3)
34      integer dlo(3), dhi(3)
35      integer g_fca, g_fcb, g_sda, g_sdb, ii
36      integer g_d1(3),g_rhs,g_u(2)
37      integer i, j, i_total
38      double precision tol2e
39      character*256 cphf_rhs, cphf_sol
40c
41      logical  cphf2, file_write_ga, file_read_ga, cphf
42      external cphf2, file_write_ga, file_read_ga, cphf
43c
44      logical     oskel
45      double precision valuea, valueb
46      data tol2e   /1.0d-16/
47c
48      double precision pifac, froth
49c
50      integer ilist(3,3)
51      data ilist /1,4,5, 4,2,6, 5,6,3/
52c
53      parameter(froth=4.0d0/3.0d0)
54      parameter(pifac=froth*3.14159265358979323846264338327950288419d0)
55c
56      oskel = .false.
57c
58c     Integral initialization
59c
60      call int_init(rtdb,1,basis)
61      call schwarz_init(geom,basis)
62      call hnd_giao_init(basis,1)
63      call scf_get_fock_param(rtdb,tol2e)
64c
65c     Total number of responses to be calculated
66c     i_pert FC + 6 * i_pert SD + 6 * i_resp SD-FC
67c
68      i_total = 7*i_pert + 6*i_resp
69c
70c     Create U matrix of dimension (nbf,nmo,3) and zero
71c     Use ahi for dimension and ahi array for chunking/blocking
72c
73      alo(1) = nbf
74      alo(2) = -1
75      alo(3) = -1
76      ahi(1) = nbf
77      ahi(2) = nopen(1)
78      ahi(3) = i_total
79      if (.not.nga_create(MT_DBL,3,ahi,'U matrix',alo,g_u(1))) call
80     &    errquit('hnd_sdfc: nga_create failed g_u',0,GA_ERR)
81      call ga_zero(g_u(1))
82      if (.not.nga_create(MT_DBL,3,ahi,'U matrix',alo,g_u(2))) call
83     &    errquit('hnd_sdfc: nga_create failed g_u',0,GA_ERR)
84      call ga_zero(g_u(2))
85c
86c     Construction of right-hand side CPHF
87c     Create CPHF array of proper dimension : (nocc*nvirt,atom)
88c
89      if(.not.ga_create(MT_DBL,2*nopen(1)*nvirt(1),i_total,'RHS',
90     &                  -1,-1,g_rhs))
91     &   call errquit('hnd_sdfc: ga_create failed g_rhs',0,GA_ERR)
92      call ga_zero(g_rhs)
93c
94c     Get FC and SD+FC in GA
95c
96      alo(1) = nbf
97      alo(2) = -1
98      alo(3) = -1
99      ahi(1) = nbf
100      ahi(2) = nbf
101      ahi(3) = i_pert+i_resp
102c
103c     Allocate memory
104c
105      if (.not.nga_create(MT_DBL,3,ahi,'fca matrix',alo,g_fca)) call
106     &    errquit('hnd_giaox: nga_create failed g_fca',0,GA_ERR)
107      if (.not.nga_create(MT_DBL,3,ahi,'fcb matrix',alo,g_fcb)) call
108     &    errquit('hnd_giaox: nga_create failed g_fca',0,GA_ERR)
109      ahi(3) = 6*(i_pert+i_resp)
110      if (.not.nga_create(MT_DBL,3,ahi,'sda matrix',alo,g_sda)) call
111     &    errquit('hnd_giaox: nga_create failed g_sda',0,GA_ERR)
112      if (.not.nga_create(MT_DBL,3,ahi,'sda matrix',alo,g_sdb)) call
113     &    errquit('hnd_giaox: nga_create failed g_sda',0,GA_ERR)
114c
115c     Calculate integrals
116c
117      call ga_zero(g_fca)
118      call ga_zero(g_sda)
119      call int_giao_1ega(basis,basis,g_fca,'fc',coords,
120     &                   i_pert+i_resp,oskel)
121      call int_giao_1ega(basis,basis,g_sda,'sd+fc',coords,
122     &                   i_pert+i_resp,oskel)
123c
124c     Take out FC part from SD+FC (i.e. from xx, yy and zz components)
125c
126      alo(1) = 1
127      ahi(1) = nbf
128      alo(2) = 1
129      ahi(2) = nbf
130      blo(1) = 1
131      bhi(1) = nbf
132      blo(2) = 1
133      bhi(2) = nbf
134      do ixy = 1, i_pert+i_resp
135         alo(3) = 1+(ixy-1)*6
136         ahi(3) = 1+(ixy-1)*6
137         blo(3) = ixy
138         bhi(3) = ixy
139         do i = 1, 3
140            call nga_add_patch(1.0d0,g_sda,alo,ahi,pifac,g_fca,blo,bhi,
141     &                         g_sda,alo,ahi)
142            alo(3) = alo(3)+1
143            ahi(3) = ahi(3)+1
144         enddo
145      enddo
146c
147c     NGA dimension arrays for copying will be the same every time
148c     Also third NGA dimension for any of the three dimensional
149c     arrays will be the same everytime (running from 1 to 3)
150c     So, lets define them once and for all in blo and bhi
151c
152      blo(1) = 1
153      bhi(1) = nopen(1)*nvirt(1)
154      blo(2) = 1
155      bhi(2) = i_pert
156c
157c     Now that all the work is done on the integrals in Alpha, copy to Beta
158c
159      call ga_copy(g_fca,g_fcb)
160      call ga_copy(g_sda,g_sdb)
161c
162c     ga_rhs(a,i) = ga_rhs(a,i) + FC(a,i)
163c     Transform FC and SD to MO and add to g_rhs
164c
165      call giao_aotomo(g_fca,vectors(1),nopen(1),nvirt(1),1,i_pert,nbf)
166      call giao_aotomo(g_fcb,vectors(2),nopen(2),nvirt(2),1,i_pert,nbf)
167      call giao_aotomo(g_sda,vectors(1),nopen(1),nvirt(1),1,
168     &                 6*(i_pert+i_resp),nbf)
169      call giao_aotomo(g_sdb,vectors(2),nopen(2),nvirt(2),1,
170     &                 6*(i_pert+i_resp),nbf)
171c
172c     Add to g_rhs
173c
174      alo(1) = nopen(1)+1
175      ahi(1) = nmo
176      alo(2) = 1
177      ahi(2) = nopen(1)
178      alo(3) = 1
179      ahi(3) = i_pert
180      call nga_add_patch(1.0d0,g_rhs,blo,bhi,0.5d0,g_fca,alo,ahi,
181     &                   g_rhs,blo,bhi)
182      blo(1) = blo(1)+nopen(1)*nvirt(1)
183      bhi(1) = bhi(1)+nopen(2)*nvirt(2)
184      call nga_add_patch(1.0d0,g_rhs,blo,bhi,-0.5d0,g_fcb,alo,ahi,
185     &                   g_rhs,blo,bhi)
186c
187      blo(1) = 1
188      bhi(1) = nopen(1)*nvirt(1)
189      blo(2) = i_pert+1
190      bhi(2) = i_pert+6*(i_pert+i_resp)
191      ahi(3) = 6*(i_pert+i_resp)
192      call nga_add_patch(1.0d0,g_rhs,blo,bhi,0.5d0,g_sda,alo,ahi,
193     &                   g_rhs,blo,bhi)
194      blo(1) = blo(1)+nopen(1)*nvirt(1)
195      bhi(1) = bhi(1)+nopen(2)*nvirt(2)
196      call nga_add_patch(1.0d0,g_rhs,blo,bhi,-0.5d0,g_sdb,alo,ahi,
197     &                   g_rhs,blo,bhi)
198      call ga_scale(g_rhs,-4.0d0)
199c
200c     Some memory cleanup here
201c
202      if (.not.ga_destroy(g_fcb)) call
203     &    errquit('hnd_sdfc: ga_destroy failed g_fcb',0,GA_ERR)
204      if (.not.ga_destroy(g_sdb)) call
205     &    errquit('hnd_sdfc: ga_destroy failed g_sdb',0,GA_ERR)
206c
207c     Write ga_rhs to disk
208c
209      call cphf_fname('cphf_rhs',cphf_rhs)
210      call cphf_fname('cphf_sol',cphf_sol)
211      if(.not.file_write_ga(cphf_rhs,g_rhs)) call errquit
212     $  ('hnd_sdfc: could not write cphf_rhs',0, DISK_ERR)
213c
214      call schwarz_tidy()
215      call int_terminate()
216c
217c     Call the CPHF routine
218c
219      if (.not.cphf2(rtdb)) call errquit
220     $  ('hnd_sdfc: failure in cphf ',0, RTDB_ERR)
221c
222c     Occ-virt blocks are the solution pieces of the CPHF
223c     Read solution vector from disk and put solutions in U matrices
224c
225      call int_init(rtdb,1,basis)
226      call schwarz_init(geom,basis)
227      call ga_zero(g_rhs)
228      if(.not.file_read_ga(cphf_sol,g_rhs)) call errquit
229     $  ('hnd_sdfc: could not read cphf_rhs',0, DISK_ERR)
230      alo(1) = nopen(1)+1
231      ahi(1) = nmo
232      alo(2) = 1
233      ahi(2) = nopen(1)
234      alo(3) = 1
235      ahi(3) = i_total
236      blo(1) = 1
237      bhi(1) = nopen(1)*nvirt(1)
238      blo(2) = 1
239      bhi(2) = i_total
240      call nga_copy_patch('n',g_rhs,blo,bhi,g_u(1),alo,ahi)
241      blo(1) = blo(1)+nopen(1)*nvirt(1)
242      bhi(1) = bhi(1)+nopen(2)*nvirt(2)
243      call nga_copy_patch('n',g_rhs,blo,bhi,g_u(2),alo,ahi)
244c
245      if (.not.ga_destroy(g_rhs)) call
246     &    errquit('hnd_sdfc: ga_destroy failed g_rhs',0,GA_ERR)
247c
248c     From U matrices, generate the perturbed density matrices D1
249c     C1 = C0 * U10
250c     D1 = [(C1*C0+) + (C0*C1+)]
251c
252      alo(1) = nbf
253      alo(2) = -1
254      alo(3) = -1
255      ahi(1) = nbf
256      ahi(2) = nbf
257      ahi(3) = i_total
258      if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1(1))) call
259     &    errquit('hnd_sdfc: nga_create failed g_d1(1)',0,GA_ERR)
260      if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1(2))) call
261     &    errquit('hnd_sdfc: nga_create failed g_d1(2)',0,GA_ERR)
262c
263      alo(1) = 1
264      alo(2) = 1
265      blo(1) = 1
266      blo(2) = 1
267      clo(1) = 1
268      chi(1) = nbf
269      clo(2) = 1
270      chi(2) = nbf
271      dlo(1) = 1
272      dhi(1) = nbf
273      dlo(2) = 1
274      do ii = 1, 2
275         dhi(2) = nopen(ii)
276      do ifld = 1, i_total
277         alo(3) = ifld
278         ahi(3) = ifld
279         blo(3) = ifld
280         bhi(3) = ifld
281         clo(3) = ifld
282         chi(3) = ifld
283         dlo(3) = ifld
284         dhi(3) = ifld
285         bhi(1) = nbf
286         bhi(2) = nmo
287         ahi(1) = nmo
288         ahi(2) = nopen(ii)
289c
290c     Make C1
291c
292         call nga_matmul_patch('n','n',1.0d0,0.0d0,vectors(ii),blo,bhi,
293     &                         g_u(ii),alo,ahi,g_d1(ii),dlo,dhi)
294         call nga_copy_patch('n',g_d1(ii),dlo,dhi,g_u(ii),dlo,dhi)
295         bhi(1) = nbf
296         bhi(2) = nopen(ii)
297         ahi(1) = nopen(ii)
298         ahi(2) = nbf
299c
300c     Make D1
301c
302         call nga_matmul_patch('n','t',1.0d0,0.0d0,vectors(ii),blo,bhi,
303     &                      g_u(ii),alo,ahi,g_d1(ii),clo,chi)
304         call nga_matmul_patch('n','t',1.0d0,0.0d0,g_u(ii),blo,bhi,
305     &                      vectors(ii),alo,ahi,g_d1(ii),clo,chi)
306      enddo
307      enddo
308c
309      call ga_sync()
310c
311      if (.not.ga_destroy(g_u(1))) call
312     &    errquit('hnd_sdfc: ga_destroy failed g_u(1)',0,GA_ERR)
313      if (.not.ga_destroy(g_u(2))) call
314     &    errquit('hnd_sdfc: ga_destroy failed g_u(2)',0,GA_ERR)
315c
316c     Now we have in g_d1(nmo,nmo,3) the derivative densities and
317c     hence we can calculate the contributions to the FC term
318c
319      call ga_zero(g_fca)
320      call int_giao_1ega(basis,basis,g_fca,'fc',coords(1,i_pert+1),
321     &                   i_resp,oskel)
322      alo(1) = 1
323      ahi(1) = nbf
324      alo(2) = 1
325      ahi(2) = nbf
326      blo(1) = 1
327      bhi(1) = nbf
328      blo(2) = 1
329      bhi(2) = nbf
330      do i = 1, ipairs
331c        alo(3) = unique(translate(i))
332c        ahi(3) = unique(translate(i))
333c        blo(3) = unique(translate(i+ipairs))
334c        bhi(3) = unique(translate(i+ipairs))
335         alo(3) = translate(i)
336         ahi(3) = translate(i)
337         blo(3) = translate(i+ipairs)
338         bhi(3) = translate(i+ipairs)
339         valuea=nga_ddot_patch(g_d1(1),'n',alo,ahi,g_fca,'n',blo,bhi)
340         valueb=nga_ddot_patch(g_d1(2),'n',alo,ahi,g_fca,'n',blo,bhi)
341c              if (ga_nodeid().eq.0)
342c    &            write(6,'(A,2I4,4F12.6)') 'FC tens i j com',
343c    &            pairlist(ixy),pairlist(ixy+ipairs),
344c    &            tensor(1,1,1,ixy),(valuea-valueb)*0.5d0*nfc,
345c    &             valuea,valueb
346         do ixy = 1, 3
347            tensor(ixy,ixy,1,i)=(valuea-valueb)*nfc*0.5d0
348         enddo
349      enddo
350c
351c     Calculate the SD-SD tensor. First calculate the sd+fc integrals
352c     and subtract the fc contribution
353c
354      call ga_zero(g_sda)
355      call int_giao_1ega(basis,basis,g_sda,'sd+fc',coords(1,i_pert+1),
356     &                   i_resp,oskel)
357      alo(1) = 1
358      alo(2) = 1
359      blo(1) = 1
360      bhi(1) = nbf
361      blo(2) = 1
362      bhi(2) = nbf
363      clo(1) = 1
364      chi(1) = nbf
365      clo(2) = 1
366      chi(2) = nbf
367      do i = 1, i_resp
368         alo(3) = 1+(i-1)*6
369         ahi(3) = 1+(i-1)*6
370         blo(3) = i
371         bhi(3) = i
372         do ixy = 1, 3
373            call nga_add_patch(1.0d0,g_sda,alo,ahi,pifac,g_fca,blo,bhi,
374     &                         g_sda,alo,ahi)
375            alo(3) = alo(3)+1
376            ahi(3) = ahi(3)+1
377         enddo
378      enddo
379c
380c     Calculate the tensor contributions to the SD-SD term
381c
382      do ixy = 1, ipairs
383c        i = unique(translate(ixy))
384c        j = unique(translate(ixy+ipairs))
385         i = translate(ixy)
386         j = translate(ixy+ipairs)
387c        RAW contributions for reference
388c        do ix = 1, 3
389c           do iy = 1, 3
390c              alo(3) = i_pert + (i-1)*6 + ilist(ix,iy)
391c              ahi(3) = i_pert + (i-1)*6 + ilist(ix,iy)
392c              blo(3) = (j-1)*6 + ilist(ix,iy)
393c              bhi(3) = (j-1)*6 + ilist(ix,iy)
394c              valuea = nga_ddot_patch(g_sda,'n',blo,bhi,
395c    &                                 g_d1(1),'n',alo,ahi)
396c              valueb = nga_ddot_patch(g_sda,'n',blo,bhi,
397c    &                                 g_d1(2),'n',alo,ahi)
398c              if (ga_nodeid().eq.0)
399c    &            write(6,'(A,4I4,5F12.6)') 'SD raw i j com',
400c    &            pairlist(ixy),pairlist(ixy+ipairs),ix,iy,
401c    &            (valuea-valueb)*nsd*0.5d0,
402c    &             valuea,valueb,nsd
403c           enddo
404c        enddo
405         do ix = 1, 3
406            do iy = 1, 3
407               valuea = 0.0d0
408               valueb = 0.0d0
409               do iz = 1, 3
410                  alo(3) = i_pert + (i-1)*6 + ilist(iz,iy)
411                  ahi(3) = i_pert + (i-1)*6 + ilist(iz,iy)
412                  blo(3) = (j-1)*6 + ilist(iz,ix)
413                  bhi(3) = (j-1)*6 + ilist(iz,ix)
414                  valuea = valuea + nga_ddot_patch(g_sda,'n',blo,bhi,
415     &                                             g_d1(1),'n',alo,ahi)
416                  valueb = valueb + nga_ddot_patch(g_sda,'n',blo,bhi,
417     &                                             g_d1(2),'n',alo,ahi)
418               enddo
419               tensor(iy,ix,2,ixy)=(valuea-valueb)*nsd*0.5d0
420            enddo
421         enddo
422      enddo
423c
424c     The final term to be calculated is the SD-FC tensorterm , which consists
425c     of two contributions:
426c     D1(SD)_N * FC_M + D1(SD)_M * FC_N  N=perturbing, M=responding
427c     We have all the contributions for the first term in the matrices,
428c     hence lets do the ddots
429c
430      do ixy = 1, ipairs
431c        i = unique(translate(ixy))
432c        j = unique(translate(ixy+ipairs))
433         i = translate(ixy)
434         j = translate(ixy+ipairs)
435         do ix = 1, 3
436            do iy = 1, 3
437               alo(3) = i_pert + (i-1)*6 + ilist(ix,iy)
438               ahi(3) = i_pert + (i-1)*6 + ilist(ix,iy)
439               blo(3) = j
440               bhi(3) = j
441               valuea = nga_ddot_patch(g_fca,'n',blo,bhi,
442     &                                 g_d1(1),'n',alo,ahi)
443               valueb = nga_ddot_patch(g_fca,'n',blo,bhi,
444     &                                 g_d1(2),'n',alo,ahi)
445               tensor(ix,iy,3,ixy)=(valuea-valueb)*0.5d0
446c              if (ga_nodeid().eq.0)
447c    &            write(6,'(A,4I4,3F12.6)') 'SDFC NM tens i j com',
448c    &            pairlist(ixy),pairlist(ixy+ipairs),
449c    &            ix,iy,tensor(ix,iy,3,ixy)*nsdfc,valuea,valueb
450            enddo
451         enddo
452c        if (ga_nodeid().eq.0) print*,''
453      enddo
454c
455c     For the D1(SD)_M * FC_N  N=perturbing, M=responding we need to calculate the
456c     fc integrals for the perturbing atoms. We already calculated the D1(SD) for
457c     the responding atoms as extra terms.
458c
459      call ga_zero(g_fca)
460      call int_giao_1ega(basis,basis,g_fca,'fc',coords,
461     &                   i_pert,oskel)
462      do ixy = 1, ipairs
463c        i = unique(translate(ixy))
464c        j = unique(translate(ixy+ipairs))
465         i = translate(ixy)
466         j = translate(ixy+ipairs)
467         do ix = 1, 3
468            do iy = 1, 3
469               valuea = 0.0d0
470               valueb = 0.0d0
471               alo(3) = 7*i_pert + (j-1)*6 + ilist(ix,iy)
472               ahi(3) = 7*i_pert + (j-1)*6 + ilist(ix,iy)
473               blo(3) = i
474               bhi(3) = i
475               valuea = valuea + nga_ddot_patch(g_fca,'n',blo,bhi,
476     &                                          g_d1(1),'n',alo,ahi)
477               valueb = valueb + nga_ddot_patch(g_fca,'n',blo,bhi,
478     &                                          g_d1(2),'n',alo,ahi)
479               tensor(ix,iy,3,ixy)=tensor(ix,iy,3,ixy)+
480     &                             (valuea-valueb)*0.5d0
481               tensor(ix,iy,3,ixy)=tensor(ix,iy,3,ixy)*nsdfc
482c              if (ga_nodeid().eq.0)
483c    &            write(6,'(A,4I4,4F12.6)') 'SDFC MN tens i j com',
484c    &            pairlist(ixy),pairlist(ixy+ipairs),
485c    &            ix,iy,tensor(ix,iy,3,ixy),(valuea-valueb)*0.5d0*nsdfc,
486c    &             valuea,valueb
487            enddo
488         enddo
489c        if (ga_nodeid().eq.0) print*,''
490      enddo
491c
492      call ga_sync()
493c
494c     Clean up all remaining memory
495c
496      if (.not.ga_destroy(g_d1(1))) call
497     &    errquit('hnd_sdfc: ga_destroy failed g_d1(1)',0,GA_ERR)
498      if (.not.ga_destroy(g_d1(2))) call
499     &    errquit('hnd_sdfc: ga_destroy failed g_d1(2)',0,GA_ERR)
500      if (.not.ga_destroy(g_sda)) call
501     &    errquit('hnd_sdfc: ga_destroy failed g_sda',0,GA_ERR)
502      if (.not.ga_destroy(g_fca)) call
503     &    errquit('hnd_sdfc: ga_destroy failed g_fca',0,GA_ERR)
504c
505      call schwarz_tidy()
506      call int_terminate()
507c
508      return
509      end
510