← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/property/giao_b1_movecs_tools.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/property/giao_b1_movecs_tools.F
 
1
      subroutine update_rhs_fock2e(
 
2
     &                      g_rhs,  ! in/out: 
 
3
     &                      vectors,! in : MO vectors
 
4
     &                      rtdb,   ! in : rtdb  handle
 
5
     &                      basis,  ! in : basis handle
 
6
     &                      geom,   ! in : geom  handle
 
7
     &                      g_dens, ! in : e-density
 
8
     &                      nocc,   ! in : nr. occ  shells
 
9
     &                      nvirt,  ! in : nr. virt shells
 
10
     &                      npol,   ! in : nr. polarizations
 
11
     &                      nbf,    ! in : nr. basis functions
 
12
     &                      nmo,    ! in : nr. MOs   
 
13
     &                      xfac,   ! in : exchange factor
 
14
     &                      tol2e,  ! in : tolerance coeff.
 
15
     &                      debug)  ! in : logical for debugging
 
16
c
 
17
c Author : Fredy W. Aquino
 
18
c Date   : 03-15-12
 
19
c Note.- Modified from original aoresponse source code
 
20
c        for extension to spin-unrestricted case
 
21
c        original aoresponse source code was written by 
 
22
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
23
c --> Experimental (not published yet)
 
24
 
 
25
       implicit none
 
26
#include "errquit.fh"
 
27
#include "global.fh"
 
28
#include "mafdecls.fh"
 
29
#include "msgids.fh"
 
30
#include "stdio.fh"
 
31
#include "geom.fh"
 
32
#include "rtdb.fh"
 
33
#include "bas.fh"
 
34
#include "prop.fh"
 
35
      integer g_fock,g_rhs
 
36
      integer rtdb,basis,geom
 
37
      integer npol,ispin
 
38
      integer vectors(npol),g_dens(3)
 
39
      integer ifld,ndir,nbf,nmo,disp,shift,      
 
40
     &        nocc(npol),nvirt(npol)
 
41
      double precision tol2e,xfac
 
42
      integer alo(3), ahi(3),
 
43
     &        blo(3), bhi(3) 
 
44
      logical debug
 
45
      external new_giao_2e,giao_aotomo
 
46
      ndir=3 ! = nr directions (x,y,z)
 
47
c     Remaining term is Perturbed (GIAO) two-electron term times
 
48
c     Unperturbed density Calculate Sum(r,s) D0(r,s) * G10(m,n,r,s) in
 
49
c     AO basis
 
50
      alo(1) = -1 
 
51
      alo(2) = -1
 
52
      alo(3) =  1
 
53
      ahi(1) = nbf
 
54
      ahi(2) = nbf
 
55
      ahi(3) = ndir*npol
 
56
      if (.not.nga_create(MT_DBL,ndir,ahi,'Fock matrix',
 
57
     &    alo,g_fock)) call 
 
58
     &    errquit('giao_b1: nga_create failed g_fock',0,GA_ERR)
 
59
      call ga_zero(g_fock)
 
60
      if(use_theory.eq.'dft') then
 
61
         ifld = 4
 
62
         if (.not. rtdb_put(rtdb,'fock_xc:calc_type',mt_int,1,ifld))
 
63
     $      call errquit('giao_b1: rtdb_put failed',0,RTDB_ERR)
 
64
      endif
 
65
 
 
66
       call new_giao_2e(geom,basis,nbf,tol2e,
 
67
     &                  g_dens, !  in: e-denstiy 
 
68
     &                  g_fock, ! out: fock matrix
 
69
     &                  xfac,
 
70
     &                  npol)
 
71
 
 
72
      if(use_theory.eq.'dft') then
 
73
         ifld = 0
 
74
         if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, ifld))
 
75
     $      call errquit('giao_b1: rtdb_put failed',0,RTDB_ERR)
 
76
         if(.not. rtdb_put(rtdb,'bgj:xc_active', MT_LOG, 1, .false.))
 
77
     $       call errquit('giao_b1: rtdb_put of xc_active failed',0,
 
78
     &       RTDB_ERR)
 
79
      endif
 
80
c
 
81
c     Transform to MO basis and add to right-hand-side
 
82
      call giao_aotomo(g_fock,vectors,nocc,nvirt,npol,ndir,nbf)
 
83
      do ispin=1,npol
 
84
       disp=ndir*(ispin-1)
 
85
       alo(1) = nocc(ispin)+1
 
86
       ahi(1) = nmo
 
87
       alo(2) = 1
 
88
       ahi(2) = nocc(ispin)
 
89
       alo(3) = disp+1
 
90
       ahi(3) = disp+ndir
 
91
       shift=nocc(1)*nvirt(1)*(ispin-1)
 
92
       blo(1) = shift+1
 
93
       bhi(1) = shift+nocc(ispin)*nvirt(ispin)
 
94
       blo(2) = 1
 
95
       bhi(2) = ndir
 
96
       call nga_add_patch(1.0d0, g_rhs,blo,bhi,
 
97
     &                    1.0d0,g_fock,alo,ahi,
 
98
     &                           g_rhs,blo,bhi)
 
99
      enddo ! end-loop-ispin
 
100
      if (debug) then
 
101
       if (ga_nodeid().eq.0)
 
102
     &  write(*,*) '------- g_fock2e-nw ---- START'
 
103
        call ga_print(g_fock)
 
104
       if (ga_nodeid().eq.0)
 
105
     &  write(*,*) '------- g_fock2e-nw ---- END'
 
106
      endif ! end-if-debug
 
107
      if (.not.ga_destroy(g_fock)) call 
 
108
     &    errquit('giao_b1: ga_destroy failed g_fock',0,GA_ERR)
 
109
      return
 
110
      end
 
111
 
 
112
      subroutine update_rhs_shfock(g_rhs,  ! in/out: RHS used for cphf2/3
 
113
     &                             g_d1,   ! in    :
 
114
     &                             vectors,! in    : MO vectors
 
115
     &                             rtdb,   ! in    : rtdb  handle
 
116
     &                             geom,   ! in    : geom  handle
 
117
     &                             basis,  ! in    : basis handle 
 
118
     &                             jfac,   ! in    : exch factors
 
119
     &                             kfac,   ! in    : exch factors
 
120
     &                             tol2e,  ! in    : tolerance coeff
 
121
     &                             nocc,   ! in    : nr occ shells
 
122
     &                             nvirt,  ! in    : nr vir shells
 
123
     &                             npol,   ! in    : nr. polarizations
 
124
     &                             nbf,    ! in    : nr. basis functions
 
125
     &                             nmo,    ! in    : nr. MOs
 
126
     &                             debug)  ! in    : =.true. for debugging
 
127
c
 
128
c Purpose: Updating g_rhs with g_fock from shell_fock_build()
 
129
c Author : Fredy W. Aquino
 
130
c Date   : 03-15-12
 
131
c Note.- Modified from original aoresponse source code
 
132
c        for extension to spin-unrestricted case
 
133
c        original aoresponse source code was written by 
 
134
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
135
c         
 
136
c --> Experimental (not published yet)
 
137
 
 
138
       implicit none
 
139
#include "errquit.fh"
 
140
#include "global.fh"
 
141
#include "mafdecls.fh"
 
142
#include "msgids.fh"
 
143
#include "stdio.fh"
 
144
#include "geom.fh"
 
145
#include "rtdb.fh"
 
146
#include "bas.fh"
 
147
#include "prop.fh"
 
148
#include "bgj.fh"
 
149
#include "case.fh"
 
150
      integer g_rhs,
 
151
     &        g_d1,g_d2,g_fock,g_s10
 
152
      integer geom,basis,rtdb
 
153
      integer npol,disp,ndens,ispin
 
154
      integer vectors(npol),nocc(npol),nvirt(npol)
 
155
      integer ifld,ndir,nbf,nmo,shift
 
156
      integer alo(3), ahi(3), 
 
157
     &        blo(3), bhi(3),
 
158
     &        clo(3), chi(3)
 
159
      logical debug
 
160
      double precision jfac(12),kfac(12),tol2e
 
161
      external shell_fock_build,
 
162
     &         shell_fock_build_cam,
 
163
     &         add_fock ! located in hnd_shift_zora.F
 
164
      ndir=3 ! nr directions (x,y,z)
 
165
c -------- Creating g_d2 ---------------START
 
166
c Note.- g_d2 =(g_d1 g_d1)
 
167
      ndens=ndir*npol
 
168
      clo(1) = ndens*2
 
169
      clo(2) = nbf
 
170
      clo(3) = nbf
 
171
      chi(1) =  1  
 
172
      chi(2) = -1 
 
173
      chi(3) = -1
 
174
      if (.not.nga_create(MT_DBL,3,clo,'g_d2 matrix',
 
175
     &                    chi,g_d2)) 
 
176
     &  call errquit('gprelim_fock: nga_create failed g_d2',
 
177
     &                0,GA_ERR)
 
178
       call ga_zero(g_d2)
 
179
       blo(1) = 1
 
180
       bhi(1) = ndens   
 
181
       blo(2) = 1
 
182
       bhi(2) = nbf
 
183
       blo(3) = 1
 
184
       bhi(3) = nbf
 
185
      do ispin=1,npol 
 
186
       disp=ndens*(ispin-1) 
 
187
       alo(1) = disp+1
 
188
       ahi(1) = disp+ndens   
 
189
       alo(2) = 1
 
190
       ahi(2) = nbf
 
191
       alo(3) = 1
 
192
       ahi(3) = nbf
 
193
       call nga_copy_patch('n',g_d1,blo,bhi,
 
194
     &                         g_d2,alo,ahi) 
 
195
      enddo ! end-loop-ispin
 
196
c -------- Creating g_d2 ---------------END
 
197
c     Build "first order fock matrix"
 
198
      if (use_theory.eq.'dft') then
 
199
         if(.not. rtdb_put(rtdb,'bgj:xc_active', MT_LOG, 1, .true.))
 
200
     $     call errquit('hess_cphf: rtdb_put of xc_active failed',0,
 
201
     &       RTDB_ERR)
 
202
         if(.not. rtdb_put(rtdb,'fock_xc:calc_type', MT_INT, 1, 2))
 
203
     $     call errquit('hess_cphf: rtdb_put of calc_type failed',0,
 
204
     &       RTDB_ERR)
 
205
         if(.not. rtdb_put(rtdb,'fock_j:derfit', MT_LOG, 1, .false.))
 
206
     $     call errquit('hess_cphf: rtdb_put of j_derfit failed',0,
 
207
     &       RTDB_ERR)
 
208
      endif
 
209
      clo(1) = ndir*npol*2
 
210
      clo(2) = nbf
 
211
      clo(3) = nbf
 
212
      chi(1) =  1  
 
213
      chi(2) = -1 
 
214
      chi(3) = -1
 
215
      if (.not.nga_create(MT_DBL,3,clo,'Fock matrix',chi,g_fock)) call 
 
216
     &    errquit('giao_b1: nga_create failed g_fock',0,GA_ERR)
 
217
      call ga_zero(g_fock)
 
218
c      if (ga_nodeid().eq.0)
 
219
c     &  write(*,*) 'cam_exch=',cam_exch
 
220
c from hnd_giaox.F
 
221
c Note: Just the exchange: jfac = 0.d0 (see above)
 
222
      if (.not.cam_exch) then
 
223
 
 
224
         call shell_fock_build(geom, basis,0,ndir*npol*2,
 
225
     $                         jfac,kfac,tol2e,
 
226
     &                         g_d2,  ! input
 
227
     &                         g_fock,! output
 
228
     &                         .false.)
 
229
 
 
230
      else
 
231
 
 
232
         call shell_fock_build_cam(geom, basis,0,ndir*npol*2,
 
233
     $                         jfac,kfac,tol2e,
 
234
     &                         g_d2,  ! input
 
235
     &                         g_fock,! output
 
236
     &                         .false.)
 
237
 
 
238
      end if
 
239
 
 
240
c      if (ga_nodeid().eq.0) then
 
241
c        write(*,70) nocc(1)   ,nocc(2),
 
242
c     &              nvirt(1)  ,nvirt(2),npol,nbf,nmo
 
243
c 70    format('BEF_add_fock:: nocc =(',i5,',',i5,') ',
 
244
c     &        'nvirt=(',i5,',',i5,') ',
 
245
c     &        '(npol,nbf,nmo)=(',i3,',',i3,',',i3,')')
 
246
c      endif
 
247
 
 
248
 
 
249
      if (debug) then
 
250
       if (ga_nodeid().eq.0)
 
251
     &  write(*,*) '------- g_fock-nw ---- START'
 
252
       call ga_print(g_fock)
 
253
       if (ga_nodeid().eq.0)
 
254
     &  write(*,*) '------- g_fock-nw ---- END'
 
255
      endif ! end-if-debug
 
256
      if(use_theory.eq.'dft') then
 
257
         if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, 0))
 
258
     $      call errquit('giaox: rtdb_put failed',0,RTDB_ERR)
 
259
      endif
 
260
c Note.- add_fock() is defined in hnd_gshift_zora.F
 
261
 
 
262
      call add_fock(g_rhs, ! out: accumulated rhs expression
 
263
     &             g_fock, !  in: Fock-term
 
264
     &            vectors, !  in: MO  coeffs
 
265
     &                nbf, !  in: nr. basis functions
 
266
     &                nmo, !  in: nr. MOs (occ+virt)
 
267
     &               npol, !  in: nr. of polarizations
 
268
     &               nocc, !  in: nr. occ     MOs
 
269
     &              nvirt) !  in: nr. virtual MOs
 
270
 
 
271
      if (debug) then
 
272
       if (ga_nodeid().eq.0)
 
273
     &  write(*,*) '---- g_rhs-AFT-shfock-inside-- START'
 
274
        call ga_print(g_rhs)
 
275
      if (ga_nodeid().eq.0)
 
276
     &  write(*,*) '---- g_rhs-AFT-shfock-inside---  END'
 
277
      endif ! end-if-debug
 
278
       if (.not.ga_destroy(g_d2)) call 
 
279
     &    errquit('giao_b1: ga_destroy failed g_fock',0,GA_ERR)
 
280
       if (.not.ga_destroy(g_fock)) call 
 
281
     &    errquit('giao_b1: ga_destroy failed g_fock',0,GA_ERR)
 
282
      return
 
283
      end
 
284
 
 
285
      subroutine get_d1_giao_b1(g_d1,   ! out:
 
286
     &                          g_u,    ! in :
 
287
     &                          vectors,! in : MO vectors
 
288
     &                          nocc,   ! in : nr. occ shells
 
289
     &                          npol,   ! in : nr. polarizations
 
290
     &                          nbf,    ! in : nr. basis functions
 
291
     &                          nmo,    ! in : nr. MOs
 
292
     &                          debug)  ! in : =.true. for debugging
 
293
c
 
294
c Purpose: get_d1 in giao_b1_movecs()
 
295
c Author : Fredy W. Aquino
 
296
c Date   : 03-15-12
 
297
c Note.- Modified from original aoresponse source code
 
298
c        for extension to spin-unrestricted case
 
299
c        original aoresponse source code was written by 
 
300
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
301
c --> Experimental (not published yet)
 
302
 
 
303
       implicit none
 
304
#include "errquit.fh"
 
305
#include "global.fh"
 
306
#include "mafdecls.fh"
 
307
#include "msgids.fh"
 
308
#include "stdio.fh"
 
309
#include "prop.fh"
 
310
      integer g_d1,g_s10
 
311
      integer ispin,npol
 
312
      integer vectors(npol),nocc(npol),
 
313
     &        g_u(npol)
 
314
      integer ifld,ndir,nbf,nmo,disp
 
315
      integer alo(3), ahi(3), 
 
316
     &        blo(3), bhi(3),
 
317
     &        clo(3), chi(3),
 
318
     &        dlo(3), dhi(3)
 
319
      logical debug
 
320
      ndir=3 ! nr directions (x,y,z)
 
321
      alo(1) = nbf
 
322
      alo(2) = -1
 
323
      alo(3) = -1
 
324
      ahi(1) = nbf
 
325
      ahi(2) = nbf
 
326
      ahi(3) = 1
 
327
      if (.not.nga_create(MT_DBL,3,ahi,'s10 matrix',alo,g_s10)) call 
 
328
     &    errquit('giao_b1: nga_create failed g_s01',0,GA_ERR)
 
329
      clo(1) = ndir*npol
 
330
      clo(2) = nbf
 
331
      clo(3) = nbf
 
332
      chi(1) =  1  
 
333
      chi(2) = -1 
 
334
      chi(3) = -1
 
335
      if (.not.nga_create(MT_DBL,3,clo,'D10 matrix',chi,g_d1)) call 
 
336
     &    errquit('giao_b1: nga_create failed g_d1',0,GA_ERR)
 
337
      call ga_zero(g_d1)
 
338
      do ispin=1,npol
 
339
       call ga_zero(g_s10)
 
340
       blo(1) = 1
 
341
       bhi(1) = nbf
 
342
       blo(2) = 1
 
343
       clo(2) = 1
 
344
       chi(2) = nbf
 
345
       clo(3) = 1
 
346
       chi(3) = nbf
 
347
       dlo(1) = 1
 
348
       dhi(1) = nbf
 
349
       dlo(2) = 1
 
350
       dhi(2) = nocc(ispin)
 
351
       disp=ndir*(ispin-1) ! for (clo,chi)
 
352
c     Create "perturbed density matrix" for closed-closed g_u block
 
353
       do ifld=1,ndir
 
354
        alo(1) = 1
 
355
        ahi(1) = nmo
 
356
        alo(2) = 1
 
357
        ahi(2) = nocc(ispin)
 
358
        alo(3) = ifld
 
359
        ahi(3) = ifld
 
360
        dlo(3) = 1
 
361
        dhi(3) = 1
 
362
        bhi(2) = nmo 
 
363
        if (debug) then
 
364
         if (ga_nodeid().eq.0) then
 
365
          write(*,17) ifld,
 
366
     &               alo(1),ahi(1),alo(2),ahi(2),
 
367
     &               alo(3),ahi(3),
 
368
     &               blo(1),bhi(1),blo(2),bhi(2),
 
369
     &               blo(3),bhi(3),
 
370
     &               dlo(1),dhi(1),dlo(2),dhi(2),
 
371
     &               dlo(3),dhi(3)
 
372
 17       format('1(',i3,')::alo-ahi=(',i3,',',i3,',',
 
373
     &          i3,',',i3,',',i3,',',i3,') ',
 
374
     &          'blo-bhi=(',i3,',',i3,',',
 
375
     &          i3,',',i3,',',i3,',',i3,') ',
 
376
     &          'dlo-dhi=(',i3,',',i3,',',
 
377
     &          i3,',',i3,',',i3,',',i3,') ')
 
378
         endif
 
379
        endif ! end-if-debug
 
380
        call nga_matmul_patch('n','n',1.0d0,0.0d0,
 
381
     &                vectors(ispin),blo,bhi,  
 
382
     &                g_u(ispin)    ,alo,ahi,
 
383
     &                g_s10         ,dlo,dhi)  
 
384
        alo(1) = 1
 
385
        ahi(1) = nocc(ispin)
 
386
        alo(2) = 1
 
387
        ahi(2) = nbf ! nmo fix lindep 05-02-12
 
388
        alo(3) = 1
 
389
        ahi(3) = 1
 
390
        bhi(2) = nocc(ispin)
 
391
        clo(1) = disp+ifld
 
392
        chi(1) = disp+ifld
 
393
c     Minus sign as we subtract it from the RHS as we do not include 
 
394
c     it in the LHS
 
395
        if (debug) then
 
396
         if (ga_nodeid().eq.0) then
 
397
          write(*,16) ifld,
 
398
     &               alo(1),ahi(1),alo(2),ahi(2),
 
399
     &               alo(3),ahi(3),
 
400
     &               blo(1),bhi(1),blo(2),bhi(2),
 
401
     &               blo(3),bhi(3),
 
402
     &               clo(1),chi(1),clo(2),chi(2),
 
403
     &               clo(3),chi(3)
 
404
 16       format('2(',i3,')::alo-ahi=(',i3,',',i3,',',
 
405
     &          i3,',',i3,',',i3,',',i3,') ',
 
406
     &          'blo-bhi=(',i3,',',i3,',',
 
407
     &          i3,',',i3,',',i3,',',i3,') ',
 
408
     &          'clo-chi=(',i3,',',i3,',',
 
409
     &          i3,',',i3,',',i3,',',i3,') ')
 
410
         endif
 
411
        endif ! end-if-debug
 
412
        call nga_matmul_patch('n','t',-1.0d0,0.0d0,
 
413
     &             vectors(ispin),blo,bhi,
 
414
     &                      g_s10,alo,ahi,
 
415
     &                       g_d1,clo,chi)  ! nbf x nbf  for (dir,ispin)
 
416
       enddo ! end-loop-ifld
 
417
      enddo ! end-loop-ispin
 
418
      if (.not.ga_destroy(g_s10)) call 
 
419
     &    errquit('giao_b1: ga_destroy failed g_s10',0,GA_ERR)
 
420
      return
 
421
      end
 
422
 
 
423
      subroutine update_rhs_threeAOints(
 
424
     &                       g_rhs,  ! in/out: RHS used for cphf2/3
 
425
     &                       vectors,! in    : MO vectors
 
426
     &                       rtdb,   ! in    : rtdb  handle 
 
427
     &                       basis,  ! in    : basis handle 
 
428
     &                       nocc,   ! in    : nr occ  shells
 
429
     &                       nvirt,  ! in    : nr virt shells
 
430
     &                       npol,   ! in    : nr. polarizations
 
431
     &                       nbf,    ! in    : nr. basis functions
 
432
     &                       nmo,    ! in    : nr. MOs
 
433
     &                       debug)  ! in    : logical for debugging
 
434
c
 
435
c Author : Fredy W. Aquino
 
436
c Date   : 03-15-12
 
437
c Note.- Modified from original aoresponse source code
 
438
c        for extension to spin-unrestricted case
 
439
c        original aoresponse source code was written by 
 
440
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
441
c --> Experimental (not published yet)
 
442
 
 
443
       implicit none
 
444
#include "errquit.fh"
 
445
#include "global.fh"
 
446
#include "mafdecls.fh"
 
447
#include "msgids.fh"
 
448
#include "stdio.fh"
 
449
#include "geom.fh"
 
450
#include "prop.fh"
 
451
#include "bas.fh"
 
452
#include "rtdb.fh"
 
453
#include "apiP.fh"
 
454
#include "bgj.fh"
 
455
#include "case.fh"
 
456
      integer npol
 
457
      integer g_rhs,g_s10,g_s10_1,vectors(npol) 
 
458
      integer rtdb,basis
 
459
      integer nocc(npol),nvirt(npol),disp,shift,ispin
 
460
      integer ifld,ndir,nbf,nmo
 
461
      integer alo(3),ahi(3), 
 
462
     &        blo(3),bhi(3)
 
463
      integer nbq,nextbq,ncosbq ! for COSMO (adding solvent effects)
 
464
      logical oskel,debug
 
465
      integer nat
 
466
      parameter (nat=1)
 
467
      double precision origin(3)
 
468
      data origin/0d0,0d0,0d0/
 
469
      external giao_aotomo
 
470
c      external geom_extbq_ncenter
 
471
c     Current CPHF does not handle symmetry 
 
472
c     Making C1 geometry and store it on rtdb
 
473
      oskel = .false.
 
474
      ndir=3 ! nr directions (x,y,z)
 
475
c     Get H10 in GA, reusing g_s10 array
 
476
      alo(1) = nbf
 
477
      alo(2) = -1
 
478
      alo(3) = -1
 
479
      ahi(1) = nbf
 
480
      ahi(2) = nbf
 
481
      ahi(3) = ndir
 
482
       if (.not.nga_create(MT_DBL,3,ahi,'s10 matrix',
 
483
     &                    alo,g_s10_1)) 
 
484
     &  call errquit('gprelim_fock: nga_create failed g_s10_1',
 
485
     &               0,GA_ERR)
 
486
       call ga_zero(g_s10_1)
 
487
      ahi(3) = ndir*npol
 
488
       if (.not.nga_create(MT_DBL,3,ahi,'s10 matrix',
 
489
     &                    alo,g_s10)) 
 
490
     &  call errquit('gprelim_fock: nga_create failed g_s10',
 
491
     &               0,GA_ERR)
 
492
      call ga_zero(g_s10)
 
493
      call int_giao_1ega(basis,basis,
 
494
     &                   g_s10_1,'l10',  ! out: g_s10 
 
495
     &                   origin,nat,oskel)
 
496
      call int_giao_1ega(basis,basis,
 
497
     &                   g_s10_1,'tv10', ! out: g_s10 updated 
 
498
     &                   origin,nat,oskel)
 
499
c
 
500
c     Get external and cosmo bq contribution
 
501
      nbq    = 0
 
502
      nextbq = 0
 
503
      ncosbq = 0
 
504
      if(geom_extbq_on()) nextbq = geom_extbq_ncenter()
 
505
      nbq = nextbq ! external bq's
 
506
      if (rtdb_get(rtdb,'cosmo:nefc',mt_int,1,ncosbq))
 
507
     &    nbq = ncosbq ! cosmo bq's
 
508
      if (nextbq.gt.0.and.ncosbq.gt.0)
 
509
     &    nbq = nextbq + ncosbq  ! tally up cosmo and external bqs
 
510
c
 
511
c     if (ga_nodeid().eq.0) write(6,*) "nbq: ", nbq
 
512
      if (nbq.gt.0) then
 
513
        call int_giao_1ega(basis,basis,
 
514
     &                     g_s10_1,'bq10', ! out
 
515
     &                     origin,nat,oskel)
 
516
      end if
 
517
c --------- g_s10_1 --> g_s10 --------- START
 
518
       blo(1) = 1
 
519
       bhi(1) = nbf ! nmo fix lindep 05-02-12
 
520
       blo(2) = 1
 
521
       bhi(2) = nbf ! nmo fix lindep 05-02-12
 
522
       blo(3) = 1
 
523
       bhi(3) = ndir
 
524
      do ispin=1,npol  
 
525
       disp=ndir*(ispin-1) 
 
526
       alo(1) = 1
 
527
       ahi(1) = nbf ! nmo fix lindep 05-02-12
 
528
       alo(2) = 1
 
529
       ahi(2) = nbf ! nmo fix lindep 05-02-12
 
530
       alo(3) = disp+1
 
531
       ahi(3) = disp+ndir  
 
532
       call nga_copy_patch('n',g_s10_1,blo,bhi,
 
533
     &                           g_s10,alo,ahi) 
 
534
      enddo ! end-loop-ispin
 
535
c --------- g_s10_1 --> g_s10 --------- END
 
536
      if (.not.ga_destroy(g_s10_1)) call 
 
537
     &    errquit('giao_b1: ga_destroy failed g_s10',0,GA_ERR)
 
538
c
 
539
c     ga_rhs(a,i) = ga_rhs(a,i) + H10(a,i)
 
540
c     Transform H10 to MO and add to g_rhs
 
541
      call giao_aotomo(g_s10,vectors,nocc,nvirt,npol,ndir,nbf)
 
542
      do ispin=1,npol
 
543
       disp=ndir*(ispin-1) 
 
544
       alo(1) = nocc(ispin)+1
 
545
       ahi(1) = nmo
 
546
       alo(2) = 1
 
547
       ahi(2) = nocc(ispin)
 
548
       alo(3) = disp+1
 
549
       ahi(3) = disp+ndir
 
550
       shift=nocc(1)*nvirt(1)*(ispin-1)
 
551
       blo(1) = shift+1
 
552
       bhi(1) = shift+nocc(ispin)*nvirt(ispin)
 
553
       blo(2) = 1
 
554
       bhi(2) = ndir
 
555
       call nga_add_patch(1.0d0,g_rhs,blo,bhi,
 
556
     &                    1.0d0,g_s10,alo,ahi,
 
557
     &                          g_rhs,blo,bhi)
 
558
      enddo ! end-loop-ispin
 
559
c      if (ga_nodeid().eq.0)
 
560
c     & write(*,*) '------- g_s10-nw ---- START'
 
561
c      call ga_print(g_s10)
 
562
c      if (ga_nodeid().eq.0)
 
563
c     & write(*,*) '------- g_s10-nw ---- END'
 
564
c
 
565
c     Cleanup g_s10 as we do not need it right now
 
566
      if (.not.ga_destroy(g_s10)) call 
 
567
     &    errquit('giao_b1: ga_destroy failed g_s10',0,GA_ERR)
 
568
      return
 
569
      end
 
570
 
 
571
      subroutine update_rhs_eS10(
 
572
     &               g_rhs,  !in/out:
 
573
     &               g_u,    !out:
 
574
     &               g_sket1,!out:
 
575
     &               eval,   !in : energy values
 
576
     &               vectors,!in : MO vectors
 
577
     &               nocc,   !in : nr.   occupied MOs
 
578
     &               nvirt,  !in : nr. unoccupied MOs
 
579
     &               npol,   !in : nr. polarizations
 
580
     &               nbf,    !in : nr. basis functions
 
581
     &               nmo,    !in : nr. MOs
 
582
     &               basis,  !in : basis handle
 
583
     &               debug)  !in : logical var for debugging
 
584
c
 
585
c Author : Fredy W. Aquino
 
586
c Date   : 03-15-12
 
587
c Note.- Modified from original aoresponse source code
 
588
c        for extension to spin-unrestricted case
 
589
c        original aoresponse source code was written by 
 
590
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
591
c --> Experimental (not published yet)
 
592
 
 
593
      implicit none
 
594
#include "errquit.fh"
 
595
#include "global.fh"
 
596
#include "mafdecls.fh"
 
597
#include "msgids.fh"
 
598
#include "stdio.fh"
 
599
#include "geom.fh"
 
600
#include "prop.fh"
 
601
#include "bas.fh"
 
602
#include "rtdb.fh"
 
603
#include "apiP.fh"
 
604
#include "bgj.fh"
 
605
#include "case.fh"
 
606
      integer basis
 
607
      integer npol,nbf,nmo     
 
608
      integer vectors(npol)
 
609
      integer nocc(npol),nvirt(npol),
 
610
     &        iocc,disp,disp1,ispin,
 
611
     &        ndir,shift,
 
612
     &        alo(3),ahi(3),
 
613
     &        blo(3),bhi(3),
 
614
     &        clo(3),chi(3)
 
615
      double precision eval(nmo*npol),toscl
 
616
      integer g_rhs,g_u(npol),g_s10,g_s10_1,g_sket1
 
617
      double precision origin(3)
 
618
      data origin/0d0,0d0,0d0/
 
619
      logical oskel,debug
 
620
      integer nat
 
621
      parameter (nat=1)
 
622
c     Current CPHF does not handle symmetry 
 
623
c     Making C1 geometry and store it on rtdb
 
624
      oskel=.false.
 
625
      ndir=3 ! nr. directions (x,y,z)
 
626
c     NGA dimension arrays for copying will be the same every time
 
627
c     Also third NGA dimension for any of the three dimensional
 
628
c     arrays will be the same everytime (running from 1 to 3)
 
629
c     So, lets define them once and for all in blo and bhi    
 
630
c     Get S10 in GA and transform to MO set (virt,occ)
 
631
      alo(1) = nbf
 
632
      alo(2) = -1
 
633
      alo(3) = -1
 
634
      ahi(1) = nbf
 
635
      ahi(2) = nbf
 
636
      ahi(3) = ndir
 
637
      if (.not.nga_create(MT_DBL,3,ahi,'s10 matrix',alo,g_s10_1)) call 
 
638
     &    errquit('giao_b1: nga_create failed g_s01',0,GA_ERR)
 
639
      call ga_zero(g_s10_1)
 
640
      ahi(3) = ndir*npol
 
641
       if (.not.nga_create(MT_DBL,3,ahi,'s10 matrix',
 
642
     &                    alo,g_s10)) 
 
643
     &  call errquit('gprelim_fock: nga_create failed g_s10',
 
644
     &               0,GA_ERR)
 
645
      call ga_zero(g_s10)
 
646
      call int_giao_1ega(basis,basis,
 
647
     &                   g_s10_1,'s10', ! out: g_s10 FA-LBL-B
 
648
     &                   origin,nat,oskel)
 
649
c -------- create g_s10 --------------- START
 
650
       blo(1) = 1
 
651
       bhi(1) = nbf ! nmo fix-lindep 05-02-12
 
652
       blo(2) = 1
 
653
       bhi(2) = nbf ! nmo fix-lindep 05-02-12
 
654
       blo(3) = 1
 
655
       bhi(3) = ndir
 
656
      do ispin=1,npol  
 
657
       disp=ndir*(ispin-1) 
 
658
       alo(1) = 1
 
659
       ahi(1) = nbf ! nmo fix-lindep 05-02-12
 
660
       alo(2) = 1
 
661
       ahi(2) = nbf ! nmo fix-lindep 05-02-12
 
662
       alo(3) = disp+1
 
663
       ahi(3) = disp+ndir
 
664
       call nga_copy_patch('n',g_s10_1,blo,bhi,
 
665
     &                         g_s10  ,alo,ahi) 
 
666
      enddo ! end-loop-ispin
 
667
       if (.not.ga_destroy(g_s10_1)) call 
 
668
     &   errquit('giao_b1: ga_destroy failed vectors',0,GA_ERR)
 
669
c -------- create g_s10 --------------- END
 
670
c After giao_aotomo valid dim(g_s10): nmo x nmo
 
671
      call giao_aotomo(g_s10,vectors,nocc,nvirt,npol,ndir,nbf) 
 
672
      if (debug) write (luout,*) 'S10 done'
 
673
c     while we are calculating integrals, let's also determine
 
674
c     the 'half' overlap derivative, used later in the calling
 
675
c     routine
 
676
      call ga_zero(g_sket1)
 
677
      call int_giao_1ega(basis,basis,
 
678
     &                   g_sket1,'srxRb', ! out: g_sket1 FA-LBL-A
 
679
     &                   origin,nat,oskel)
 
680
      if (debug) write (luout,*) 'S1-ket done'
 
681
c     g_sket1 will not be used further here. It is one of the 
 
682
c     output results of this routine.
 
683
c     Broceed with the computation of the B-field perturbed
 
684
c     MO coefficients.
 
685
      do ispin=1,npol
 
686
c     ga_rhs(a,i) = ga_rhs(a,i) - e(i) * S10(a,i)
 
687
c     Scale (occ,virt) block g_s10 with - (minus) eigenvalues 
 
688
       disp1=ndir*(ispin-1)
 
689
       alo(1) = nocc(ispin)+1
 
690
       ahi(1) = nmo
 
691
       alo(3) = disp1+1
 
692
       ahi(3) = disp1+ndir
 
693
       clo(1) = 1
 
694
       chi(1) = nocc(ispin)
 
695
       clo(2) = 1
 
696
       chi(2) = nocc(ispin)
 
697
       clo(3) = 1
 
698
       chi(3) = ndir
 
699
c       disp   = nmo*(ispin-1)
 
700
       disp   = nbf*(ispin-1) ! fix-lindep 05-02-12
 
701
       do iocc=1,nocc(ispin)
 
702
        alo(2) = iocc
 
703
        ahi(2) = iocc
 
704
        toscl  =-eval(disp+iocc) 
 
705
        call nga_scale_patch(g_s10,alo,ahi,toscl) 
 
706
       enddo ! end-loop-iocc
 
707
c     Copy to ga_rhs 
 
708
c     alo(1) and ahi(1) the same as before
 
709
       alo(2) = 1
 
710
       ahi(2) = nocc(ispin)
 
711
       shift=nocc(1)*nvirt(1)*(ispin-1)
 
712
       blo(1) = shift+1
 
713
       bhi(1) = shift+nocc(ispin)*nvirt(ispin)
 
714
       blo(2) = 1
 
715
       bhi(2) = ndir
 
716
       call nga_copy_patch('n',g_s10,alo,ahi,
 
717
     &                         g_rhs,blo,bhi)
 
718
      
 
719
c     Construct occ-occ part of the three U matrices
 
720
c     Occ-occ blocks for each field direction are defined as -1/2 S10
 
721
c     Scale (occ,occ) block g_s10 with -1/2 and add to g_u
 
722
c     alo(2) and ahi(2) will stay as 1 and nclosed(1) for a while
 
723
       alo(1) = 1
 
724
       ahi(1) = nocc(ispin)
 
725
       call nga_scale_patch(g_s10,alo,ahi,-0.5d0)
 
726
       call nga_copy_patch('n',g_s10     ,alo,ahi,
 
727
     &                         g_u(ispin),clo,chi)
 
728
      enddo ! end-loop-ispin
 
729
       if (debug) write (luout,*) 'S10 in occ-occ done'
 
730
       if (.not.ga_destroy(g_s10)) call 
 
731
     &   errquit('giao_b1: ga_destroy failed vectors',0,GA_ERR)
 
732
      return
 
733
      end