← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/ddscf/rohf_hessv3_ext.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/ddscf/rohf_hessv3_ext.F
 
1
c ... jochen: this comes from rohf_hessv2.F but has everything
 
2
c     "v2" replaced by "v3", for use with frequency-dependent response
 
3
c     and is called from the cphf_solve3 part of the cphf code
 
4
c
 
5
c ... jochen: further mods were made to accomodate the situation that
 
6
c     we might have damping in the response. That causes all quantities to
 
7
c     have an imaginary part, too
 
8
 
 
9
      subroutine rohf_hessv3_ext(acc, 
 
10
     &                       g_x, g_ax, 
 
11
     &                       g_x_im, g_Ax_im,
 
12
     &   omega, limag, lifetime, gamwidth, ncomp)
 
13
c
 
14
c     Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
15
c     Date   : 03-24-12
 
16
c     Note.- Equivalent to rohf_hessv3() (routine written by 
 
17
c            J. Autschbach) but using optimizing routines.
 
18
 
 
19
      implicit none
 
20
#include "errquit.fh"
 
21
#include "crohf.fh"
 
22
#include "cscf.fh"
 
23
#include "stdio.fh"
 
24
#include "util.fh"
 
25
#include "global.fh"
 
26
c     
 
27
c     $Id: rohf_hessv3.F 19983 2011-02-21 02:53:21Z niri $
 
28
c
 
29
c ... jochen: these two arrays now have two components:
 
30
      integer g_x(2)  ! [input]  A-matrix elements for density matrix
 
31
      integer g_ax(2) ! [output] Perturbed Fock operator
 
32
c ... jochen: also, we might have imaginary components:
 
33
      integer g_x_im(2)  ! [input]  A-matrix elements, Im
 
34
      integer g_ax_im(2) ! [output] Perturbed Fock operator, Im
 
35
 
 
36
      double precision acc, omega, gamwidth
 
37
      logical limag, lifetime
 
38
c     
 
39
      integer gtype,grow,gcol,growp,gcolp, ipm, ncomp
 
40
      logical oprint, debug
 
41
      external rohf_hessv_xx3_ext
 
42
c
 
43
c     ================================================================
 
44
 
 
45
      debug = (.false. .and. ga_nodeid().eq.0) ! for code development
 
46
 
 
47
      if (debug) write (6,*)
 
48
     &   'rohf_hessv3: limag, omega, lifetime, gamwidth',
 
49
     &   limag, omega, lifetime, gamwidth
 
50
c
 
51
c     Check for debug 
 
52
c     
 
53
      oprint= util_print('rohf_hessv2',print_debug)
 
54
      if (crohf_init_flag.ne.1)
 
55
     $     call errquit('rohf_hessv3: ROHF internal block invalid',0,
 
56
     &       UNKNOWN_ERR)
 
57
c
 
58
c ... jochen: use first component for the dimension checks.
 
59
c     the second component MUST have the same dimensions
 
60
c     otherwise there will be problems
 
61
      call ga_inquire(g_x(1),gtype,grow,gcol)
 
62
      if (grow.ne.crohf_vlen)
 
63
     $     call errquit('rohf_hessv3: invalid vector length',0,
 
64
     &       UNKNOWN_ERR)
 
65
      call ga_inquire(g_ax(1),gtype,growp,gcolp)
 
66
      if (growp.ne.crohf_vlen)
 
67
     $     call errquit('rohf_hessv3: invalid vector length',0,
 
68
     &       UNKNOWN_ERR)
 
69
      if (gcol.ne.gcolp)
 
70
     $     call errquit('rohf_hessv3: invalid no. of vectors',0,
 
71
     &       UNKNOWN_ERR)
 
72
c     
 
73
c     Call internal routine
 
74
c  
 
75
      if (debug) write (6,*) 'calling rohf_hessv_xx3'
 
76
 
 
77
      call rohf_hessv_xx3_ext( basis, geom, nbf, nmo,
 
78
     $     nclosed, nopen,
 
79
     $     pflg, 
 
80
     &     g_movecs, 
 
81
     &     oskel, noskew,
 
82
     $     crohf_g_fcv, crohf_g_fpv, crohf_g_fcp,
 
83
     $     acc, lshift, 
 
84
     &     g_x, g_ax, 
 
85
     &     g_x_im, g_Ax_im, 
 
86
     &     omega, limag,
 
87
     &     lifetime, gamwidth, ncomp)
 
88
 
 
89
      if (debug) write (6,*) 'back from rohf_hessv_xx3'
 
90
c
 
91
c     Zap numbers much smaller than acc to ensure hard zeroes 
 
92
c     remain unpolluted ... cannot use a threshold larger than the
 
93
c     integral accuracy since can break symmetry in non-abelian groups
 
94
c     Also must ensure that the threshold tends to zero to permit
 
95
c     tight convergence.
 
96
c
 
97
c ... jochen: screen components
 
98
      do ipm = 1,ncomp
 
99
        call ga_screen(g_ax(ipm),
 
100
     &       max(min(acc*acc,acc*0.01d0,1d-12),1d-16))
 
101
        if (lifetime) call ga_screen(g_ax_im(ipm), 
 
102
     &       max(min(acc*acc,acc*0.01d0,1d-12),1d-16))
 
103
      enddo
 
104
c
 
105
      end
 
106
 
 
107
      subroutine rohf_hessv_xx3_ext( basis, geom, 
 
108
     &     nbf, nmo, nclosed, nopen, 
 
109
     $     pflg,
 
110
     $     g_movecs, oskel, noskew, g_fcv, g_fpv, g_fcp,
 
111
     $     acc, lshift, 
 
112
     &     g_x, g_ax, 
 
113
     &     g_x_im, g_Ax_im, 
 
114
     &     omega, limag,
 
115
     &     lifetime, gamwidth, ncomp)
 
116
c     Note.- Improvements to this subroutine done by
 
117
c            Fredy W. Aquino, Northwestern University (Oct 2012)
 
118
c            Using rohf_hessv_2e2_opt() instead of old rohf_hessv_2e2()
 
119
c                  rohf_hessv_2e3_opt() instead of old rohf_hessv_2e3() 
 
120
 
 
121
C     $Id: rohf_hessv3.F 19983 2011-02-21 02:53:21Z niri $
 
122
      implicit none
 
123
#include "errquit.fh"
 
124
#include "global.fh"
 
125
#include "mafdecls.fh"
 
126
#include "rtdb.fh"
 
127
#include "bgj.fh"
 
128
c     
 
129
      integer basis, geom
 
130
      integer nbf, nmo, nclosed, nopen
 
131
      integer pflg
 
132
      integer g_movecs
 
133
      logical oskel, noskew
 
134
      integer g_fcv, g_fpv, g_fcp
 
135
      double precision acc
 
136
      double precision lshift
 
137
c ... jochen: input arrays g_x and g_Ax have two components here
 
138
      integer g_x(2), g_ax(2), vlen, nvec, g_tmp, gtype, ipm
 
139
      integer g_x_im(2)  ! [input]  A-matrix elements, Im
 
140
      integer g_ax_im(2) ! [output] Perturbed Fock operator, Im
 
141
 
 
142
      double precision omega, gamwidth, wls, wlsim,omg(2)
 
143
      logical limag, lifetime
 
144
      integer ncomp
 
145
      logical debug
 
146
      external rohf_hessv_2e3_opt 
 
147
 
 
148
c
 
149
c     =================================================================
 
150
 
 
151
      debug = (.false. .and. ga_nodeid().eq.0) ! for code development
 
152
 
 
153
c      debug = (.true. .and. ga_nodeid().eq.0) ! for code development
 
154
c
 
155
      if (debug) write (6,*) 'hessv3: omega =',omega
 
156
      if (debug) write (6,*) 'hessv3: limag =',limag
 
157
      if (debug) write (6,*)
 
158
     &   'hessv3: lifetime, gamwidth =',lifetime, gamwidth
 
159
c
 
160
      do ipm = 1,ncomp
 
161
        call ga_zero(g_Ax(ipm))
 
162
        if (lifetime) call ga_zero(g_Ax_im(ipm))
 
163
      end do
 
164
      if (pflg.gt.2 .or. pflg.le.0) then
 
165
         call errquit('rohf_hessv_xx: pflg invalid ', pflg,
 
166
     &       UNKNOWN_ERR)
 
167
      endif
 
168
c
 
169
c ... jochen: to be consistent with the preconditioner, where
 
170
c     the level shift is added, we need to do the same thing here
 
171
c     and also add and subtract the frequency times 4 (it is times
 
172
c     4 because of the factors of 4 in rohf_hessv_1e and in the
 
173
c     preconditioner)
 
174
c     During a response calculation, pflg is equal to 2
 
175
c
 
176
c     what do we do here? Compare Gauss' paper Eqs. (32) and (135):
 
177
c     The lhs of the CPHF equations contain a term
 
178
c     (e_a - e_i -/+ omega) U_ai. First, we initialize g_Ax with
 
179
c     the term proportional to omega, then we add the delta-e term
 
180
c     (the e's are the orbital energies, calculated in hessv_1e as
 
181
c     the diagonal of the Fock matrix transformed to the MO basis)
 
182
 
 
183
      if (pflg .gt. 0) then
 
184
        omg(1)=-omega
 
185
        omg(2)= omega
 
186
        if (.not.lifetime) then
 
187
c        no damping: initialize Ax with terms proportional omega
 
188
         do ipm=1,ncomp
 
189
          wls = lshift + 4d0 * omg(ipm)
 
190
          call ga_dadd(wls,g_x(ipm),0.d0,g_ax(ipm),g_ax(ipm))
 
191
         enddo ! end-loop-ipm
 
192
        else                    ! lifetime
 
193
c        take care of damping here: Re and Im are coupled by gamwidth
 
194
         do ipm=1,ncomp
 
195
          wls   = lshift + 4d0* omg(ipm)
 
196
          wlsim = -4d0 * gamwidth
 
197
c          write(*,19) lshift,wls,wlsim
 
198
c 19       format('RE:(lshift,wls,wlsim)=(',
 
199
c     &           f15.8,',',f15.8,',',f15.8,')') 
 
200
 
 
201
          call ga_dadd(wls  ,g_x(ipm),
 
202
     &                 wlsim,g_x_im(ipm),
 
203
     &                       g_ax(ipm))
 
204
          wls   =  4d0 * gamwidth
 
205
          wlsim = lshift + 4d0* omg(ipm)
 
206
 
 
207
c          write(*,20) lshift,wls,wlsim
 
208
c 20       format('IM:(lshift,wls,wlsim)=(',
 
209
c     &           f15.8,',',f15.8,',',f15.8,')') 
 
210
 
 
211
          call ga_dadd(wls  ,g_x(ipm),
 
212
     &                 wlsim,g_x_im(ipm),
 
213
     &                       g_ax_im(ipm))  
 
214
         enddo ! end-lopp-ipm
 
215
        endif                   ! .not.lifetime
 
216
        call ga_sync()
 
217
        if (debug) write (6,*) 'calling rohf_hessv_1e'
 
218
 
 
219
c ============== debug g_ax ==================== START
 
220
        if (debug) then
 
221
          do ipm=1,2
 
222
            if (ga_nodeid().eq.0)        
 
223
     &      write(*,*) '------- g_ax_re-0-c(',ipm,')------ START' 
 
224
            call ga_print(g_ax(ipm))
 
225
            if (ga_nodeid().eq.0)
 
226
     &      write(*,*) '------- g_ax_re-0-c(',ipm,')------ END'
 
227
          enddo ! end-loop-ipm
 
228
         if (lifetime) then
 
229
          do ipm=1,2
 
230
           if (ga_nodeid().eq.0)        
 
231
     &     write(*,*) '------- g_ax_im-0-c(',ipm,')------ START' 
 
232
           call ga_print(g_ax_im(ipm))
 
233
           if (ga_nodeid().eq.0)
 
234
     &     write(*,*) '------- g_ax_im-0-c(',ipm,')------ END'
 
235
          enddo ! end-loop-ipm
 
236
         endif ! end-if-lifetime
 
237
        endif ! end-if-debug
 
238
c ============== debug g_ax ==================== END
 
239
c 
 
240
c       next: add (e_a - e_i) times A (also called U) matrix to Ax
 
241
        do ipm=1,ncomp
 
242
         call rohf_hessv_1e(basis,geom,        ! in : handles
 
243
     &                      nmo,nclosed,nopen, ! in : (nmo,nocc) nopen=0 for closed shell
 
244
     $                      g_fcv,g_fpv,g_fcp, ! in : densities
 
245
     $                      g_x(ipm),          ! in : g_x
 
246
     &                      g_ax(ipm))         ! out: 1e contrib to Ax
 
247
        enddo ! end-loop-ipm
 
248
        if (lifetime) then
 
249
        do ipm=1,ncomp
 
250
         call rohf_hessv_1e(basis,geom,        ! in : handles
 
251
     &                      nmo,nclosed,nopen, ! in : (nmo,nocc) nopen=0 for closed shell
 
252
     $                      g_fcv,g_fpv,g_fcp, ! in : densities
 
253
     $                      g_x_im(ipm),       ! in : g_x_im
 
254
     &                      g_ax_im(ipm))      ! out: 1e contrib to Ax_im
 
255
        enddo ! end-loop-ipm
 
256
        endif ! end-if-lifetime
 
257
      endif                     ! pflg.gt.0
 
258
c ============== debug g_ax ==================== START
 
259
      if (debug) then
 
260
       do ipm=1,2
 
261
         if (ga_nodeid().eq.0)        
 
262
     &    write(*,*) '------- g_ax_re-0-d(',ipm,')------ START' 
 
263
         call ga_print(g_ax(ipm))
 
264
         if (ga_nodeid().eq.0)
 
265
     &    write(*,*) '------- g_ax_re-0-d(',ipm,')------ END'
 
266
       enddo ! end-loop-ipm
 
267
       if (lifetime) then
 
268
       do ipm=1,2
 
269
         if (ga_nodeid().eq.0)        
 
270
     &    write(*,*) '------- g_ax_im-0-d(',ipm,')------ START' 
 
271
         call ga_print(g_ax_im(ipm))
 
272
         if (ga_nodeid().eq.0)
 
273
     &    write(*,*) '------- g_ax_im-0-d(',ipm,')------ END'
 
274
       enddo ! end-loop-ipm
 
275
       endif ! end-if-lifetime
 
276
      endif ! end-if-debug
 
277
c ============== debug g_ax ==================== END
 
278
      if (pflg .gt. 1) then
 
279
c
 
280
c       the next call basically uses the current guess for the solution
 
281
c       vector x (in g_x, which is the perturbed density matrix in the
 
282
c       MO basis) and calculates the perturbed Fock operator in the MO basis.
 
283
c       real and imaginary part of that Fock operator can be handled
 
284
c       separately here
 
285
 
 
286
        if (ncomp.gt.1) then    ! call 2e code for dynamic case
 
287
 
 
288
c          if (ga_nodeid().eq.0)
 
289
c     &     write(*,*) 'FA-BEF rohf_hessv_2e3_opt-RE-IM'
 
290
 
 
291
           call rohf_hessv_2e3_opt(
 
292
     &                  basis,   ! in: basis handle
 
293
     &                  geom,    ! in: geom  handle
 
294
     &                  nbf,     ! in: nr. basis functions
 
295
     &                  nmo,     ! in: nr. MOs
 
296
     &                  nclosed, ! in: nr. double occupied MOs
 
297
     &                  nopen,   ! in: nr. single occupied MOs
 
298
     $                  g_movecs,! in: MO coefficients
 
299
     &                  oskel,   ! in: =.true. ->
 
300
     &                  noskew,  ! in: =.true. -> symmetric density matrix
 
301
     &                  g_x,     ! in: 
 
302
     &                  g_x_im,  ! in: 
 
303
     &                  g_ax,    ! in/out: Hessian product
 
304
     &                  g_ax_im, ! in/out: Hessian product
 
305
     &                  acc,     ! in: accuracy Fock construction
 
306
     &                  limag,   ! in: =.true. -> imaginary component allowed
 
307
     &                  lifetime)! in: =.true. -> RE-IM =.false -> RE
 
308
 
 
309
c            if (ga_nodeid().eq.0)
 
310
c     &       write(*,*) 'FA-AFT rohf_hessv_2e3_opt-RE-IM'
 
311
 
 
312
        else                    ! call static 2e code
 
313
 
 
314
c          if (ga_nodeid().eq.0)
 
315
c     &     write(*,*) 'FA-BEF rohf_hessv_2e2_opt'
 
316
 
 
317
             call rohf_hessv_2e2_opt(
 
318
     &                              basis,     ! in : basis handle
 
319
     &                              geom,      ! in : geom  handle
 
320
     &                              nbf,       ! in : nr. basis functions
 
321
     &                              nmo,       ! in : nr. MOs vecs
 
322
     &                              nclosed,   ! in : nr. occupied MOs 
 
323
     &                              nopen,     ! in : nr. open shells (unpaired e's)
 
324
     $                              g_movecs,  ! in : MO vec coeffs 
 
325
     &                              oskel,     ! in :
 
326
     &                              noskew,    ! in : symm density ?
 
327
     &                              acc,       ! in : accuracy Fock construction   
 
328
     &                              g_x(1),    ! in : RE g_x 
 
329
     &                              g_x_im(1), ! in : IM g_x
 
330
     &                              g_ax(1),   ! in/ou : RE g_ax Hessian product
 
331
     &                              g_ax_im(1),! in/ou : IM g_ax Hessian product
 
332
     &                              lifetime)
 
333
 
 
334
c          if (ga_nodeid().eq.0)
 
335
c     &     write(*,*) 'FA-AFT rohf_hessv_2e2_opt'
 
336
 
 
337
        endif                   ! ncomp
 
338
        
 
339
      endif                     ! pflg.gt.1
 
340
c     
 
341
      end
 
342
 
 
343
      subroutine rohf_hessv3_cmplx(
 
344
     &                  acc,       ! in : accuracy
 
345
     &                  g_z,       ! in : z
 
346
     &                  g_Az1,     ! in : Az1
 
347
     &                  nsub,      ! in
 
348
     &                  omega,     ! in :
 
349
     &                  limag,     ! in :
 
350
     &                  lifetime,  ! in : 
 
351
     &                  gamwidth,  ! in :
 
352
     &                  ncomp,     ! in : nr. components (+/-)
 
353
     &                  iter_cphf) ! in: iteration nr. in cphf
 
354
c
 
355
c     Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
356
c     Date   : 03-24-12
 
357
c     Note.-  Modifying rohf_hessv3 to handle complex vars
 
358
c --> Experimental (not published yet)
 
359
 
 
360
      implicit none
 
361
#include "errquit.fh"
 
362
#include "crohf.fh"
 
363
#include "cscf.fh"
 
364
#include "stdio.fh"
 
365
#include "util.fh"
 
366
#include "global.fh"
 
367
c     
 
368
c     $Id: rohf_hessv3.F 19983 2011-02-21 02:53:21Z niri $
 
369
      integer ncomp
 
370
      integer g_z(ncomp), ! [input]  A-matrix elements for density matrix
 
371
     &        g_Az(ncomp) ! [output] Perturbed Fock operator
 
372
      integer g_Az1,
 
373
     &        nsub
 
374
      double precision acc, omega, gamwidth
 
375
      logical limag, lifetime    
 
376
      integer gtype,grow,gcol,
 
377
     &        growp,gcolp,ipm,iter_cphf
 
378
      logical oprint, debug
 
379
      external rohf_hessv_xx3_cmplx
 
380
 
 
381
      debug = (.false. .and. ga_nodeid().eq.0) ! for code development
 
382
      if (debug) write (6,*)
 
383
     &   'rohf_hessv3: limag, omega, lifetime, gamwidth',
 
384
     &   limag, omega, lifetime, gamwidth
 
385
 
 
386
c
 
387
c     Check for debug 
 
388
c     
 
389
      oprint= util_print('rohf_hessv2',print_debug)
 
390
      if (crohf_init_flag.ne.1)
 
391
     $     call errquit('rohf_hessv3: ROHF internal block invalid',0,
 
392
     &       UNKNOWN_ERR)
 
393
 
 
394
c
 
395
c ... jochen: use first component for the dimension checks.
 
396
c     the second component MUST have the same dimensions
 
397
c     otherwise there will be problems
 
398
      call ga_inquire(g_z(1),gtype,grow,gcol) 
 
399
 
 
400
c      if (ga_nodeid().eq.0) then
 
401
c       write(*,117) grow,gcol,crohf_vlen
 
402
c117    format('(grow,gcol,crohf_vlen)=(',
 
403
c     &         i4,',',i4,',',i4,')')      
 
404
c      endif
 
405
 
 
406
      if (grow.ne.crohf_vlen)
 
407
     $     call errquit('rohf_hessv3_cmplx: invalid vector length-1',
 
408
     &                  0,UNKNOWN_ERR)
 
409
 
 
410
c      call ga_inquire(g_Az(1),gtype,growp,gcolp)
 
411
 
 
412
c      if (ga_nodeid().eq.0) then
 
413
c       write(*,118) growp,gcolp,crohf_vlen
 
414
c118    format('(growp,gcolp,crohf_vlen)=(',
 
415
c     &         i4,',',i4,',',i4,')')      
 
416
c      endif
 
417
 
 
418
      goto 177 ! skip-for-the moment
 
419
      if (growp.ne.crohf_vlen)
 
420
     $     call errquit('rohf_hessv3_cmplx: invalid vector length-2',
 
421
     &                  0,UNKNOWN_ERR)
 
422
      if (gcol.ne.gcolp)
 
423
     $     call errquit('rohf_hessv3_cmplx: invalid no. of vectors',
 
424
     &                  0,UNKNOWN_ERR)
 
425
 177  continue
 
426
c     
 
427
c     Call internal routine
 
428
c  
 
429
      if (debug) write (6,*) 'calling rohf_hessv_xx3_cmplx'
 
430
 
 
431
c      if (ga_nodeid().eq.0)
 
432
c     &  write(*,*) 'BEF rohf_hessv_xx3_cmplx'
 
433
 
 
434
      call rohf_hessv_xx3_cmplx(
 
435
     &     g_z,         ! in :
 
436
     &     g_Az1,
 
437
     &     nsub,
 
438
     &     ncomp,       ! in : nr. components
 
439
     &     basis,       ! in : basis handle
 
440
     &     geom,        ! in : geom  hande
 
441
     &     nbf,         ! in : nr. basis functions
 
442
     &     nmo,         ! in : nr. MOs
 
443
     $     nclosed,     ! in : nr. occupied MOs
 
444
     &     nopen,       ! in : nr. 1/2 occupied MOs (=0 for closed shells)
 
445
     $     pflg, 
 
446
     &     g_movecs,    ! in : MO coefficients 
 
447
     &     oskel, 
 
448
     &     noskew,
 
449
     $     crohf_g_fcv, ! in : closed-virtual density
 
450
     &     crohf_g_fpv, ! in : open-virtual   density (not used)
 
451
     &     crohf_g_fcp, ! in : closed-open    density (not used)
 
452
     $     acc, 
 
453
     &     lshift, 
 
454
     &     omega,
 
455
     &     limag,
 
456
     &     lifetime,
 
457
     &     gamwidth,
 
458
     &     iter_cphf)  ! in : iteration nr. in cphf cycle
 
459
 
 
460
c      if (ga_nodeid().eq.0)
 
461
c     &  write(*,*) 'AFT rohf_hessv_xx3_cmplx'
 
462
 
 
463
      if (debug) write (6,*) 'back from rohf_hessv_xx3_cmplx'
 
464
c
 
465
c     Zap numbers much smaller than acc to ensure hard zeroes 
 
466
c     remain unpolluted ... cannot use a threshold larger than the
 
467
c     integral accuracy since can break symmetry in non-abelian groups
 
468
c     Also must ensure that the threshold tends to zero to permit
 
469
c     tight convergence.
 
470
c
 
471
c ... jochen: screen components
 
472
c ++++++++++++++++++++++++++++++++++++++++++
 
473
c =====> WARNING-UNFINISHED ROUTINE <=======
 
474
c ++++++++++++++++++++++++++++++++++++++++++
 
475
c     I need to adapt it for complex (later) --> zapping routine
 
476
c      do ipm = 1,ncomp
 
477
c        call ga_screen(g_ax(ipm),
 
478
c     &                 max(min(acc*acc,acc*0.01d0,1d-12),1d-16))
 
479
c        if (lifetime) call ga_screen(g_ax_im(ipm), 
 
480
c     &       max(min(acc*acc,acc*0.01d0,1d-12),1d-16))
 
481
c      enddo
 
482
 
 
483
      end
 
484
 
 
485
      subroutine rohf_hessv_xx3_cmplx(
 
486
     &               g_z,     ! in :
 
487
     &               g_Az1,   ! ou : product: A x z
 
488
     &               nsub,    ! in : nr. subspace
 
489
     &               ncomp,   ! in : nr. components
 
490
     &               basis,   ! in : basis handle
 
491
     &               geom,    ! in : geom  hande
 
492
     &               nbf,     ! in : nr. basis functions
 
493
     &               nmo,     ! in : nr. MOs
 
494
     &               nclosed, ! in : nr. occupied MOs
 
495
     &               nopen,   ! in : nr. 1/2 occupied MOs (=0 for closed shells)
 
496
     $               pflg,
 
497
     $               g_movecs,! in : MO coefficients 
 
498
     &               oskel, 
 
499
     &               noskew, 
 
500
     &               g_fcv,   ! in : closed-virtual density
 
501
     &               g_fpv,   ! in : open-virtual   density (not used)
 
502
     &               g_fcp,   ! in : closed-open    density (not used)
 
503
     $               acc, 
 
504
     &               lshift, 
 
505
     &               omega, 
 
506
     &               limag,
 
507
     &               lifetime, 
 
508
     &               gamwidth,
 
509
     &               iter_cphf)
 
510
c
 
511
c     Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
512
c     Date   : 03-24-12
 
513
c     -> modified from rohf_hessv_xx3()
 
514
c     $Id: rohf_hessv3.F 19983 2011-02-21 02:53:21Z niri $
 
515
c --> Experimental (not published yet)
 
516
 
 
517
      implicit none
 
518
#include "errquit.fh"
 
519
#include "global.fh"
 
520
#include "mafdecls.fh"
 
521
#include "rtdb.fh"
 
522
#include "bgj.fh"
 
523
      integer ipm,ncomp,iter_cphf
 
524
      integer g_z(ncomp), ! [input]  A-matrix elements
 
525
     &        g_Az(ncomp) ! [output] Perturbed Fock operator 
 
526
      integer g_Az1,nsub,
 
527
     &        alo(3),ahi(3),
 
528
     &        n1,maxsub,
 
529
     &        m1,m2,p1,p2,
 
530
     &        pretty
 
531
      integer g_movecs   
 
532
      integer g_fcv,
 
533
     &        g_fpv,
 
534
     &        g_fcp
 
535
      integer g_xreim,g_Axreim ! scratch GA arrays
 
536
      integer basis,geom
 
537
      integer nbf,nmo,nvir,
 
538
     &        nclosed,nopen,nreim,n
 
539
      integer pflg
 
540
      double precision acc,lshift
 
541
      integer vlen, nvec,g_tmp,gtype
 
542
      double precision omega,gamwidth,
 
543
     &                 omg(ncomp), 
 
544
     &                 wls,wlsim
 
545
      double complex wls_cmplx
 
546
      logical oskel, noskew
 
547
      logical limag,lifetime
 
548
      logical debug
 
549
      external rohf_hessv_2e3_opt_cmplx,
 
550
     &         rohf_hessv_2e2_opt_cmplx,
 
551
     &         getreorim,getreorim1,
 
552
     &         conv2complex3,conv2complex4
 
553
c
 
554
c     =================================================================
 
555
 
 
556
      debug = (.false. .and. ga_nodeid().eq.0) ! for code development
 
557
c      debug = (.true. .and. ga_nodeid().eq.0) ! for code development
 
558
 
 
559
      if (debug) write (6,*) 'hessv3: omega =',omega
 
560
      if (debug) write (6,*) 'hessv3: limag =',limag
 
561
      if (debug) write (6,*)
 
562
     &   'hessv3: lifetime, gamwidth =',lifetime, gamwidth
 
563
 
 
564
c ----- Clear sub-block ---- START
 
565
        call ga_inquire(g_Az1,gtype,n1,maxsub)
 
566
c        if (ga_nodeid().eq.0) then
 
567
c         write(*,103) n1,maxsub
 
568
c 103     format('(n1,maxsumb)=(',i4,',',i4,')')
 
569
c        endif
 
570
c
 
571
       call ga_inquire(g_z(1),gtype,n,nvec) ! get (n1,nvec)
 
572
c       write(*,178) n,nvec
 
573
c  178  format('(n,nvec)=(',i3,',',i3,')')
 
574
c
 
575
        alo(1)=1
 
576
        ahi(1)=n1
 
577
        alo(2)=nsub+1
 
578
        ahi(2)=nsub+nvec
 
579
        call nga_zero_patch(g_Az1,alo,ahi)
 
580
c ----- Clear sub-block ---- END
 
581
c
 
582
      if (pflg.gt.2 .or. pflg.le.0) then
 
583
         call errquit('rohf_hessv_xx: pflg invalid ', pflg,
 
584
     &       UNKNOWN_ERR)
 
585
      endif
 
586
c
 
587
c ... jochen: to be consistent with the preconditioner, where
 
588
c     the level shift is added, we need to do the same thing here
 
589
c     and also add and subtract the frequency times 4 (it is times
 
590
c     4 because of the factors of 4 in rohf_hessv_1e and in the
 
591
c     preconditioner)
 
592
c     During a response calculation, pflg is equal to 2
 
593
c
 
594
c     what do we do here? Compare Gauss' paper Eqs. (32) and (135):
 
595
c     The lhs of the CPHF equations contain a term
 
596
c     (e_a - e_i -/+ omega) U_ai. First, we initialize g_Ax with
 
597
c     the term proportional to omega, then we add the delta-e term
 
598
c     (the e's are the orbital energies, calculated in hessv_1e as
 
599
c     the diagonal of the Fock matrix transformed to the MO basis)
 
600
      if (pflg .gt. 0) then
 
601
        omg(1)=-omega
 
602
        omg(2)= omega
 
603
c        take care of damping here: Re and Im are coupled by gamwidth
 
604
         m1=1
 
605
         m2=n
 
606
         p1=nsub+1
 
607
         p2=nsub+nvec
 
608
c         write(*,179) m1,m2,p1,p2
 
609
c 179     format('(m1,m2,p1,p2)=(',i3,',',i3,',',i3,',',i3,')')
 
610
         do ipm=1,ncomp
 
611
          wls   = lshift + 4.0d0 * omg(ipm)
 
612
          wlsim = -4d0 * gamwidth
 
613
          wls_cmplx=dcmplx(wls,-wlsim)
 
614
c          write(*,19) lshift,wls,wlsim,wls_cmplx
 
615
c 19       format('(lshift,wls,wlsim,wls_cmplx)=(',
 
616
c     &           f15.8,',',f15.8,',',f15.8,
 
617
c     &           ',[',f15.8,',',f15.8,'])') 
 
618
c ----- copy to patch of g_Az1 --------- START
 
619
          call ga_copy_patch('n',g_z(ipm),1 ,n ,1 ,nvec,
 
620
     &                           g_Az1   ,m1,m2,p1,p2)
 
621
c ----- copy to patch of g_Az1 --------- END
 
622
c ----- Show printout of patch to be updated only --- START
 
623
c          if (ga_nodeid().eq.0)
 
624
c     &      write(*,*) '-------g_Az1(',ipm,')----START'
 
625
c            call ga_print(g_Az1)
 
626
c          if (ga_nodeid().eq.0)
 
627
c     &      write(*,*) '-------g_Az1(',ipm,')----END'
 
628
c ----- Show printout of patch to be updated only --- END
 
629
c ----- scale patch of g_Az1 --------- START
 
630
          call ga_scale_patch(g_Az1,m1,m2,p1,p2,wls_cmplx)
 
631
          m1=m1+n
 
632
          m2=m2+n
 
633
c ----- scale patch of g_Az1 --------- END
 
634
         enddo ! end-lopp-ipm
 
635
c ============== debug g_ax ==================== START
 
636
        if (debug) then
 
637
          m1=1
 
638
          m2=n
 
639
          p1=nsub+1
 
640
          p2=nsub+nvec
 
641
          do ipm=1,ncomp
 
642
            if (ga_nodeid().eq.0)        
 
643
     &      write(*,*) '------- g_Az1-c(',ipm,')------ START' 
 
644
            call ga_print(g_Az1)
 
645
            if (ga_nodeid().eq.0)
 
646
     &      write(*,*) '------- g_Az1-c(',ipm,')------ END'
 
647
           m1=m1+n
 
648
           m2=m2+n
 
649
          enddo ! end-loop-ipm
 
650
        endif ! end-if-debug
 
651
c ============== debug g_ax ==================== END
 
652
c 
 
653
c       next: add (e_a - e_i) times A (also called U) matrix to Ax
 
654
       call ga_inquire(g_z(1),gtype,n,nvec) ! get (n,nvec)
 
655
        if (.not. ga_create(MT_DBL,n,nvec,
 
656
     &     'hessv_xx3_cmplx: g_xreim',0,0,g_xreim))
 
657
     $   call errquit('hessv_xx3_cmplx: failed allocating g_xreim', 
 
658
     &                n,GA_ERR)     
 
659
        if (.not. ga_create(MT_DBL,n,nvec,
 
660
     &     'hessv_xx3_cmplx: g_xreim',0,0,g_Axreim))
 
661
     $   call errquit('hessv_xx3_cmplx: failed allocating g_xreim', 
 
662
     &                n,GA_ERR)  
 
663
 
 
664
       nvir = nmo - nclosed - nopen   
 
665
c       write(*,89) nvir,nclosed,nmo,nopen
 
666
c  89   format('(nvir,nclosed,nmo,nopen)=(',
 
667
c     &        i3,',',i3,',',i3,',',i3,')')
 
668
 
 
669
       do nreim=1,2 ! loop in RE,IM
 
670
        do ipm=1,ncomp
 
671
         call getreorim(g_xreim,  ! out : real or im arr
 
672
     &                  g_z(ipm), ! in  : = complx(g_xre,g_xim)
 
673
     &                  nvir,     ! in  : nr. virtual  MOs
 
674
     &                  nclosed,  ! in  : nr. occupied MOs
 
675
     &                  nreim)    ! in  : =1 -> re =2 -> im
 
676
 
 
677
c         if (ga_nodeid().eq.0)
 
678
c     &    write(*,*) '----BEF g_z(',ipm,',',nreim,')---START'
 
679
c        call ga_print(g_z(ipm))
 
680
c         if (ga_nodeid().eq.0)
 
681
c     &    write(*,*) '----BEF g_z(',ipm,',',nreim,')---END'
 
682
c         if (ga_nodeid().eq.0)
 
683
c     &    write(*,*) '----BEF g_xreim(',ipm,',',nreim,')---START'
 
684
c        call ga_print(g_xreim)
 
685
c         if (ga_nodeid().eq.0)
 
686
c     &    write(*,*) '----BEF g_xreim(',ipm,',',nreim,')---END'
 
687
 
 
688
         call getreorim1(g_Axreim,! out : real or im arr
 
689
     &                   g_Az1,   ! in  : = complx(g_xre,g_xim)
 
690
     &                   nsub,    ! in  : subblock index
 
691
     &                   ipm,     ! in  : = 1,2 to access slctd component
 
692
     &                   nvir,    ! in  : nr. virtual  MOs
 
693
     &                   nclosed, ! in  : nr. occupied MOs
 
694
     &                   nreim)   ! in  : =1 -> re =2 -> im
 
695
 
 
696
c         if (ga_nodeid().eq.0)
 
697
c     &    write(*,*) '----BEF g_Axreim(',ipm,',',nreim,')---START'
 
698
c         call ga_print(g_Axreim)
 
699
c         if (ga_nodeid().eq.0)
 
700
c     &    write(*,*) '----BEF g_Axreim(',ipm,',',nreim,')---END'
 
701
 
 
702
         call rohf_hessv_1e(
 
703
     &                  basis,geom,        ! in : handles
 
704
     &                  nmo,nclosed,nopen, ! in : (nmo,nocc) nopen=0 for closed shell
 
705
     $                  g_fcv,g_fpv,g_fcp, ! in : densities
 
706
     &                  g_xreim,
 
707
     &                  g_Axreim)
 
708
 
 
709
c         if (ga_nodeid().eq.0)
 
710
c     &    write(*,*) '----AFT g_Axreim(',ipm,',',nreim,')---START'
 
711
c         call ga_print(g_Axreim)
 
712
c         if (ga_nodeid().eq.0)
 
713
c     &    write(*,*) '----AFT g_Axreim(',ipm,',',nreim,')---END'
 
714
c ++++++++ update g_Az +++++++++ START
 
715
 
 
716
         call conv2complex4(g_Az1,   ! out: = history matrix complex
 
717
     &                      g_Axreim,! in : real      arr
 
718
     &                      nsub,    ! in  : subblock index
 
719
     &                      ipm,     ! in  : = 1,2 to access slctd component
 
720
     &                      nvir,    ! in  : nr. virtual  MOs
 
721
     &                      nclosed, ! in  : nr. occupied MOs
 
722
     &                      nreim)   ! in  : =1 -> re =2 -> im
 
723
 
 
724
c ++++++++ update g_Az +++++++++ END
 
725
        enddo ! end-loop-ipm
 
726
       enddo ! end-loop-nreim
 
727
        if (.not. ga_destroy(g_xreim))  call errquit
 
728
     &     ('hessv_xx3_cmplx: g_xreim',0, GA_ERR)
 
729
        if (.not. ga_destroy(g_Axreim)) call errquit
 
730
     &     ('hessv_xx3_cmplx: g_xreim',0, GA_ERR)
 
731
      endif                     ! pflg.gt.0
 
732
c ============== debug g_ax ==================== START
 
733
      if (debug) then
 
734
            if (ga_nodeid().eq.0)        
 
735
     &      write(*,*) '------- g_Az1-d------ START' 
 
736
            call ga_print(g_Az1)
 
737
            if (ga_nodeid().eq.0)
 
738
     &      write(*,*) '------- g_Az1-d------ END'
 
739
      endif ! end-if-debug
 
740
c ============== debug g_ax ==================== END
 
741
      if (pflg .gt. 1) then
 
742
c
 
743
c       the next call basically uses the current guess for the solution
 
744
c       vector x (in g_x, which is the perturbed density matrix in the
 
745
c       MO basis) and calculates the perturbed Fock operator in the MO basis.
 
746
c       real and imaginary part of that Fock operator can be handled
 
747
c       separately here
 
748
 
 
749
        if (ncomp.gt.1) then    ! call 2e code for dynamic case
 
750
 
 
751
c          if (ga_nodeid().eq.0)
 
752
c     &     write(*,*) 'FA-BEF rohf_hessv_2e3_opt_cmplx'
 
753
 
 
754
           call rohf_hessv_2e3_opt_cmplx(
 
755
     &                  g_z,      ! in :
 
756
     &                  g_Az1,    ! out: history of g_Az
 
757
     &                  nsub,     ! in :
 
758
     &                  basis,    ! in : basis handle
 
759
     &                  geom,     ! in : geom  handle
 
760
     &                  nbf,      ! in : nr. basis functions
 
761
     &                  nmo,      ! in : nr. MOs
 
762
     &                  nclosed,  ! in : nr. double occupied MOs
 
763
     &                  nopen,    ! in : nr. single occupied MOs
 
764
     $                  g_movecs, ! in : MO coefficients
 
765
     &                  oskel,    ! in : =.true. ->
 
766
     &                  noskew,   ! in : =.true. -> symmetric density matrix
 
767
     &                  acc,      ! in : accuracy Fock construction
 
768
     &                  limag,    ! in : =.true. -> imaginary component allowed
 
769
     &                  lifetime) ! in : =.true. -> RE-IM =.false -> RE
 
770
cold     &                  iter_cphf)! in: iteration nr. in cphf cycle
 
771
 
 
772
c            if (ga_nodeid().eq.0)
 
773
c     &       write(*,*) 'FA-AFT rohf_hessv_2e3_opt_cmplx'
 
774
 
 
775
c ++++++++++++++++++++++++++++++++++++++++++++
 
776
c             if (ga_nodeid().eq.0)
 
777
c     &        write(*,*) 'FA-stop-test-rohf_hessv'
 
778
c             stop
 
779
c ++++++++++++++++++++++++++++++++++++++++++++
 
780
        else                    ! call static 2e code
 
781
 
 
782
c          if (ga_nodeid().eq.0)
 
783
c     &     write(*,*) 'FA-BEF rohf_hessv_2e2_opt_cmplx'
 
784
 
 
785
           call rohf_hessv_2e2_opt_cmplx(
 
786
     &                          g_z(1),  ! in :
 
787
     &                          g_Az1,   ! out: history of g_Az
 
788
     &                          nsub,    ! in :
 
789
     &                          basis,   ! in : basis handle
 
790
     &                          geom,    ! in : geom  handle
 
791
     &                          nbf,     ! in : nr. basis functions
 
792
     &                          nmo,     ! in : nr. MOs vecs
 
793
     &                          nclosed, ! in : nr. occupied MOs 
 
794
     &                          nopen,   ! in : nr. open shells (unpaired e's)
 
795
     $                          g_movecs,! in : MO vec coeffs 
 
796
     &                          oskel,   ! in :
 
797
     &                          noskew,  ! in : symm density ?
 
798
     &                          acc,     ! in : accuracy Fock construction   
 
799
     &                          lifetime)
 
800
 
 
801
c          if (ga_nodeid().eq.0)
 
802
c     &     write(*,*) 'FA-AFT rohf_hessv_2e2_opt_cmplx'
 
803
 
 
804
        endif                   ! ncomp
 
805
        
 
806
      endif                     ! pflg.gt.1
 
807
c     
 
808
      end
 
809
 
 
810
      subroutine rohf_hessv_2e3_opt_cmplx(
 
811
     &                  g_z,
 
812
     &                  g_Az1,   ! in: (n1,maxsub) history of Az matrix (large matrix)
 
813
     &                  nsub,    ! in: point to (n1,nvec) block to be updated in g_Az1
 
814
     &                  basis,   ! in: basis handle
 
815
     &                  geom,    ! in: geom  handle
 
816
     &                  nbf,     ! in: nr. basis functions
 
817
     &                  nmo,     ! in: nr. MOs
 
818
     &                  nclosed, ! in: nr. double occupied MOs
 
819
     &                  nopen,   ! in: nr. single occupied MOs
 
820
     $                  g_movec, ! in: MO coefficients
 
821
     &                  oskel,   ! in: =.true. ->
 
822
     &                  noskew,  ! in: =.true. -> symmetric density matrix
 
823
     &                  acc,     ! in: accuracy Fock construction
 
824
     &                  limag,   ! in: =.true. -> imaginary component allowed
 
825
     &                  lifetime)! in: =.true. -> RE-IM =.false -> RE
 
826
c
 
827
c   Purpose: Optimization of rohf_hessv_2e3()
 
828
c   Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
829
c   Date   : 03-28-12
 
830
c   Note1.- Modifying rohf_hessv_2e3, reducing computation of
 
831
c          two-electron integrals by putting together 
 
832
c          symmetric+antisymmetric perturbed densities: g_dens(ipm) ipm=1,2
 
833
c          and doing a single call to shell_fock_build() when using
 
834
c          the routine shell_fock_build2().
 
835
c  Note2.- size(g_Az1)=(n1,maxsub)  n1=ncomp*n  maxsumb=maxiter*nvec
 
836
c          ncomp=2 (+/-) n=nvir*nocc maxiter=10 (usually) nvec=3 (x,y,z)
 
837
c          nsub, will point to next (n1,nvec) block to be updated 
 
838
c          --> The purpose of using g_Az1 is to reduce usage of memory.
 
839
c              by doing it I skip using g_Az(ipm) ipm=2 (n1,nvec) complex matrix.
 
840
c --> Experimental (not published yet)
 
841
 
 
842
      implicit none
 
843
#include "errquit.fh"
 
844
#include "mafdecls.fh"
 
845
#include "global.fh"
 
846
#include "util.fh"
 
847
#include "cscfps.fh"
 
848
#include "rtdb.fh"
 
849
#include "bgj.fh"
 
850
#include "stdio.fh"
 
851
#include "case.fh"
 
852
c     
 
853
c     Return the ROHF orbital 2e-Hessian vector product, g_ax = A * g_x
 
854
c
 
855
c ... jochen: modified version of rohf_hessv_2e2 which keeps track
 
856
c     of two sets of input vectors that couple via the density matrix.
 
857
c     one could likely save some memory here by re-using temp arrays
 
858
c ... jochen: Also made modifications to calculate imaginary terms due
 
859
c     to finite lifetime damping   
 
860
ccccccccccccccc This code does NOT work for open shell!!!!!ccccccccccccccccc
 
861
      integer g_z(2),g_Az(2)
 
862
      integer g_Az1,nsub
 
863
      integer m1,m2,p1,p2,pretty
 
864
      integer basis, geom       ! basis & geom handle
 
865
      integer nclosed,nvir,nopen! Basis size and occupation
 
866
      integer nbf,nmo           ! No. of linearly dependent MOs
 
867
      integer g_movec           ! MO coefficients
 
868
      logical oskel
 
869
 
 
870
      double precision acc      ! Accuracy of "Fock" construction
 
871
      logical limag             ! imaginary perturbation?    
 
872
      integer voff,xoff
 
873
      integer nmul,gtype,
 
874
     &        ivec,vlen,nvec
 
875
      integer g_tmp,g_tmp1,
 
876
     &        g_dens(4),g_fock(4),
 
877
     &        g_xreim(2)
 
878
      double precision tol2e
 
879
      logical odebug,oprint,noskew
 
880
      integer dims(3),chunk(3), 
 
881
     &        alo(3),ahi(3), 
 
882
     &        blo(2),bhi(2)
 
883
      integer ga_create_atom_blocked
 
884
      external ga_create_atom_blocked,
 
885
     &         get_undosymm_fock,update_ax_fock,
 
886
     &         shell_fock_build2,getreorim,
 
887
     &         update_gz_reorim,
 
888
     &         update_gz_reorim1 
 
889
      double precision one,zero,mone,four,half,mhalf,two,mtwo
 
890
      double precision itol_floor, itol_ceil
 
891
      parameter(itol_floor=1.d-15, itol_ceil=1.d-3)
 
892
      parameter (one=1.0d0, mone=-1.0d0, zero=0.0d0, four=4.0d0)
 
893
      parameter (half=0.5d0, mhalf=-0.5d0, two=2.0d0, mtwo=-2.0d0)
 
894
      integer ipm ! counter for density matrix components
 
895
      character*(255) cstemp
 
896
      integer g_pmats(2),g_pmata(2),g_h1mat(2),
 
897
     &        cnt,ind,indx(2,2),
 
898
     &        npol,nset,nblock,ncomp
 
899
      double precision tenm6,coef(2,2)
 
900
      parameter (tenm6 = 1d-6)
 
901
      logical lifetime,debug
 
902
      data npol /1/ ! for restricted calculations
 
903
      data indx /1,2, ! indx(1,1),indx(1,2)
 
904
     &           3,4/ ! indx(2,1),indx(2,2)
 
905
      ncomp=2 ! using two components
 
906
      call ga_inquire(g_z(1),gtype,vlen,nvec) ! out: nvec,vlen
 
907
      do ipm = 1,ncomp
 
908
       if (.not. ga_create(MT_DBL,vlen,nvec, 
 
909
     &      'hessv_2e3_opt_cmplx: g_xreim',0,0,g_xreim(ipm)))
 
910
     $   call errquit('rhessv_2e3_opt_cmplx: failed alloc g_xreim',
 
911
     &                nvec,GA_ERR)
 
912
      enddo ! end-loop-ipm
 
913
      nmul=1
 
914
      if (npol.eq. 2) nmul=2
 
915
      nset  =1 ! for RE          
 
916
      nblock=2 ! for RE
 
917
      if (lifetime) then
 
918
      nset  =2 ! for RE-IM
 
919
      nblock=4 ! for RE-IM
 
920
      endif
 
921
 
 
922
      if (oskel) 
 
923
     $   call errquit('rohf_h2e3: no way',0, UNKNOWN_ERR)
 
924
 
 
925
      debug = (.false. .and. ga_nodeid().eq.0) ! for code development
 
926
 
 
927
c      debug=.true. ! for debugging
 
928
 
 
929
      if (debug) then
 
930
       if (ga_nodeid().eq.0) then
 
931
        write(*,2001) limag,lifetime,nset,nblock
 
932
 2001   format('(limag,lifetime,nset,nblock)=(',
 
933
     &           L1,',',L1,',',i3,',',i3,')')
 
934
       endif
 
935
      endif ! end-if-debug
 
936
 
 
937
      oprint= util_print('rohf_hessv2',print_debug)
 
938
 
 
939
      if (nopen.ne.0) call errquit
 
940
     $     ('rohf_h2e3: does not work for open shells',nopen,
 
941
     &       UNKNOWN_ERR)
 
942
      odebug = util_print('rohf_hessv', print_debug)    
 
943
      tol2e = min(max(acc,itol_floor),itol_ceil)
 
944
      nvir = nmo - nclosed - nopen
 
945
      voff = nclosed + nopen + 1
 
946
 
 
947
      dims(1)  = nbf
 
948
      dims(2)  = nbf
 
949
      chunk(1) = dims(1)
 
950
      chunk(2) = -1
 
951
 
 
952
c ------ create scratch GA arrays (g_pmats,g_mata,g_h1mat---START
 
953
      do ipm = 1,ncomp
 
954
        write(cstemp,'(a,i1)') 'pmats_',ipm
 
955
        if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk,
 
956
     &     g_pmats(ipm))) call 
 
957
     &     errquit('rohf_h2e3: nga_create failed '//cstemp(1:7),
 
958
     &     0,GA_ERR)
 
959
        call ga_zero(g_pmats(ipm))
 
960
        write(cstemp,'(a,i1)') 'pmata_',ipm
 
961
 
 
962
c        if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk,
 
963
c     &     g_pmata(ipm))) call 
 
964
c     &     errquit('rohf_h2e3: nga_create failed '//cstemp(1:7),
 
965
c     &     0,GA_ERR)
 
966
c        call ga_zero(g_pmata(ipm))
 
967
 
 
968
        write(cstemp,'(a,i1)') 'h1mat_',ipm
 
969
        if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk,
 
970
     &     g_h1mat(ipm))) call 
 
971
     &     errquit('rohf_h2e3: nga_create failed '//cstemp(1:7),
 
972
     &     0,GA_ERR)
 
973
        call ga_zero(g_h1mat(ipm))
 
974
      enddo ! end-loop-ipm
 
975
c ------ create scratch GA arrays (g_pmats,g_mata,g_h1mat---END
 
976
      dims(1)  = nvec
 
977
      dims(2)  = nbf
 
978
      dims(3)  = nbf
 
979
      chunk(1) = dims(1)
 
980
      chunk(2) = -1
 
981
      chunk(3) = -1   
 
982
      call ga_sync()  
 
983
      do ipm = 1,nblock ! =2 or 4
 
984
c ... allocate g_dens=[g_dens_re g_dens_im]
 
985
        if (.not. nga_create (MT_DBL,3,dims,'CPKS dens',chunk,
 
986
     &     g_dens(ipm)))
 
987
     &     call errquit('rohf_h2e3: could not allocate g_dens',555,
 
988
     &     GA_ERR)
 
989
        call ga_zero(g_dens(ipm))
 
990
c ... allocate g_fock=[g_fock_re g_fock_im]
 
991
        if (.not. nga_create (MT_DBL,3,dims,'Fockv',chunk,
 
992
     &     g_fock(ipm)))
 
993
     &     call errquit('rohf_h2e3: could not allocate g_fock',555,
 
994
     &     GA_ERR)
 
995
        call ga_zero(g_fock(ipm))
 
996
      enddo ! end-loop-ipm
 
997
        if (debug) then
 
998
         if (ga_nodeid().eq.0)
 
999
     &    write(*,*) 'BEF get_dens_reorim-RE'
 
1000
        endif ! end-if-debug
 
1001
c ---- Copy g_z --> g_x_reim ------ START
 
1002
        do ipm=1,ncomp
 
1003
         call ga_zero(g_xreim(ipm))
 
1004
         call getreorim(g_xreim(ipm),! out : real or im arr
 
1005
     &                  g_z(ipm),    ! in  : = complx(g_xre,g_xim)
 
1006
     &                  nvir,        ! in  : nr. virtual  MOs
 
1007
     &                  nclosed,     ! in  : nr. occupied MOs
 
1008
     &                  1)           ! in  : =1 -> re =2 -> im
 
1009
        enddo ! end-loop-ipm
 
1010
c ---- Copy g_z --> g_x_reim ------ END
 
1011
      call get_dens_reorim(
 
1012
     &                    g_dens, ! in/ou: perturbed density matrix
 
1013
     &                    1,      ! in   : =1 1st block RE
 
1014
     &                    g_xreim,! in   : RE
 
1015
     &                    g_movec,! in   : MO coefficients
 
1016
     &                    nbf,    ! in   : nr. basis functions
 
1017
     &                    nmo,    ! in   : nr. MOs
 
1018
     &                    1,      ! in   : for ipol=1 -> istart=1
 
1019
     &                    nclosed,! in   : nr. occupied MOs
 
1020
     &                    nvir,   ! in   : nr. virtual  MOs 
 
1021
     &                    nvec,   ! in   : nr. directions (x,y,z)
 
1022
     &                    npol,   ! in   : nr. polarizations
 
1023
     &                    limag,  ! in   : = .true. imaginary allowed   
 
1024
     &                    g_pmats,! in   : scratch GA array
 
1025
     &                    g_pmata,! in   : scratch GA array - NOT USED 
 
1026
     &                    g_h1mat)! in   : scratch GA array   
 
1027
         if (debug) then
 
1028
          if (ga_nodeid().eq.0)
 
1029
     &     write(*,*) 'BEF get_dens_reorim-IM'
 
1030
         endif ! end-if-debug
 
1031
 
 
1032
      if (lifetime) then
 
1033
c ---- Copy g_z --> g_x_reim ------ START
 
1034
        do ipm=1,ncomp
 
1035
         call ga_zero(g_xreim(ipm))
 
1036
         call getreorim(g_xreim(ipm),! out : real or im arr
 
1037
     &                  g_z(ipm),    ! in  : = complx(g_xre,g_xim)
 
1038
     &                  nvir,        ! in  : nr. virtual  MOs
 
1039
     &                  nclosed,     ! in  : nr. occupied MOs
 
1040
     &                  2)           ! in  : =1 -> re =2 -> im
 
1041
        enddo ! end-loop-ipm
 
1042
c ---- Copy g_z --> g_x_reim ------ END
 
1043
 
 
1044
       call get_dens_reorim(
 
1045
     &                    g_dens, ! in/ou: perturbed density matrix
 
1046
     &                    2,      ! in   : =2 2nd block IM
 
1047
     &                    g_xreim,! in   : IM
 
1048
     &                    g_movec,! in   : MO coefficients
 
1049
     &                    nbf,    ! in   : nr. basis functions
 
1050
     &                    nmo,    ! in   : nr. MOs
 
1051
     &                    1,      ! in   : for ipol=1 -> istart=1
 
1052
     &                    nclosed,! in   : nr. occupied MOs
 
1053
     &                    nvir,   ! in   : nr. virtual  MOs 
 
1054
     &                    nvec,   ! in   : nr. directions (x,y,z)
 
1055
     &                    npol,   ! in   : nr. polarizations
 
1056
     &                    limag,  ! in   : = .true. imaginary allowed   
 
1057
     &                    g_pmats,! in   : scratch GA array
 
1058
     &                    g_pmata,! in   : scratch GA array - NOT USED   
 
1059
     &                    g_h1mat)! in   : scratch GA array  
 
1060
 
 
1061
      endif ! end-if-lifetime
 
1062
      do ipm = 1,ncomp
 
1063
       if (.not.ga_destroy(g_xreim(ipm))) call errquit(
 
1064
     &     'rohf_hessv3: ga_destroy failed g_xreim',0,GA_ERR)  
 
1065
       if (.not.ga_destroy(g_h1mat(ipm)))  call errquit(
 
1066
     &      'rohf_hessv3: ga_destroy failed g_h1mat',0,GA_ERR)     
 
1067
      enddo ! end-loop-ipm
 
1068
c         if (ga_nodeid().eq.0) 
 
1069
c     &     write(*,*) 'FA-BEF shell_fock_build2 ...'
 
1070
 
 
1071
      if (debug) then
 
1072
        do ipm=1,nblock
 
1073
         if (ga_nodeid().eq.0)
 
1074
     &    write(*,*) '--------g_dens(',ipm,')-------START'
 
1075
         call ga_print(g_dens(ipm))
 
1076
         if (ga_nodeid().eq.0)
 
1077
     &    write(*,*) '--------g_dens(',ipm,')-------END'
 
1078
        enddo ! end-loop-ipm
 
1079
      endif ! end-if-debug
 
1080
 
 
1081
      call shell_fock_build2(g_fock,   ! out: Fock    matrices
 
1082
     &                       g_dens,   ! in : density matrices
 
1083
     &                       geom,     ! in : geom  handle
 
1084
     &                       basis,    ! in : basis handle
 
1085
     &                       nbf,      ! in : nr. basis functions
 
1086
     &                       nvec,     ! in : nr. vecs (x,y,z)
 
1087
     &                       npol,     ! in : npol=1 for R-DFT =2 for U-DFT
 
1088
     &                       ncomp,    ! in : nr. components
 
1089
     &                       nblock,   ! in : nr. of g_dens,g_fock blocks
 
1090
     &                       .true.,   ! in : =.true. for symm dens
 
1091
     &                       tol2e,    ! in :
 
1092
     &                       debug)    ! in : = .true. -> debugging printouts
 
1093
 
 
1094
c       if (ga_nodeid().eq.0) 
 
1095
c     &  write(*,*) 'FA-AFT shell_fock_build2 ...'
 
1096
 
 
1097
      if (debug) then
 
1098
       do ipm=1,nblock
 
1099
         if (ga_nodeid().eq.0)        
 
1100
     &    write(*,*) '------- g_fock-0(',ipm,')------ START' 
 
1101
         call ga_print(g_fock(ipm))
 
1102
         if (ga_nodeid().eq.0)
 
1103
     &    write(*,*) '------- g_fock-0(',ipm,')------ END'
 
1104
       enddo ! end-loop-ipm
 
1105
      endif ! end-if-debug
 
1106
 
 
1107
      call get_undosymm_fock(
 
1108
     &            g_fock,  ! in/ou: fock matrix
 
1109
     &            nset,    ! in   : =1 g_x is real, =2 g_x is complex (g_x_re,g_x_im)
 
1110
     &            nvec,    ! in   : nr. directions (x,y,z)
 
1111
     &            nbf,     ! in   : nr. basis functions
 
1112
     &            npol,    ! in   : nr. polarizations
 
1113
     &            nmul,    ! in   : =1 npol=1 =2 npol=2 (acc. JK terms)
 
1114
     &            g_pmats, ! in   : scratch GA array
 
1115
     &            limag)   ! in   : =.true. imaginary comp. exists
 
1116
 
 
1117
c ------- Remove GA arrays:
 
1118
      do ipm = 1,ncomp
 
1119
       if (.not.ga_destroy(g_pmats(ipm))) call errquit(
 
1120
     &     'rohf_hessv3: ga_destroy failed g_pmats',0,GA_ERR)       
 
1121
c       if (.not.ga_destroy(g_pmata(ipm))) call errquit(
 
1122
c     &      'rohf_hessv3: ga_destroy failed g_pmata',0,GA_ERR)
 
1123
      enddo ! end-loop-ipm
 
1124
 
 
1125
      if (debug) then
 
1126
       do ipm=1,nblock
 
1127
         if (ga_nodeid().eq.0)        
 
1128
     &    write(*,*) '------- g_fock-unsym(',ipm,')------ START' 
 
1129
         call ga_print(g_fock(ipm))
 
1130
         if (ga_nodeid().eq.0)
 
1131
     &    write(*,*) '------- g_fock-unsym(',ipm,')------ END'
 
1132
       enddo ! end-loop-ipm
 
1133
      endif ! end-if-debug
 
1134
        
 
1135
      if (debug) then
 
1136
          if (ga_nodeid().eq.0)        
 
1137
     &    write(*,*) '------- g_Az1-0------ START' 
 
1138
          call ga_print(g_Az1)
 
1139
          if (ga_nodeid().eq.0)
 
1140
     &    write(*,*) '------- g_Az1-0------ END'
 
1141
      endif ! end-if-debug
 
1142
c     
 
1143
c     start loop over components of perturbing field
 
1144
c   
 
1145
      g_tmp = ga_create_atom_blocked(geom, basis,'rohf_h2e3: tmp')
 
1146
      g_tmp1= ga_create_atom_blocked(geom, basis,'rohf_h2e3: tmp1')
 
1147
      alo(2) = 1
 
1148
      ahi(2) = nbf
 
1149
      alo(3) = 1
 
1150
      ahi(3) = nbf
 
1151
      blo(1) = 1
 
1152
      bhi(1) = nbf
 
1153
      blo(2) = 1
 
1154
      bhi(2) = nclosed  
 
1155
      do cnt=1,nset
 
1156
       do ivec = 1, nvec
 
1157
        alo(1) = ivec
 
1158
        ahi(1) = ivec     
 
1159
        do ipm = 1,ncomp ! loop over Fock matrix components +/- here     
 
1160
          ind=indx(ipm,cnt)
 
1161
 
 
1162
          if (debug) then
 
1163
            if (ga_nodeid().eq.0) then
 
1164
             write(*,117) cnt,ivec,ipm,ind
 
1165
 117         format('XX:(cnt,ivec,ipm,ind)=(',
 
1166
     &              i3,',',i3,',',i3,',',i3,')')
 
1167
            endif
 
1168
          endif ! end-if-debug
 
1169
c         
 
1170
c          P      =  4(ij|kl) - (ik|jl) - (il|kj)
 
1171
c           ij,kl
 
1172
c     
 
1173
c          K      =  (ik|jl) + (il|kj)
 
1174
c           ij,kl
 
1175
c     
 
1176
c          cv         cv          pv   cp
 
1177
c          Z   =  2P.[D  ]  +  P.[D  + D  ]
 
1178
c     
 
1179
c          pv          cv           cp   pv
 
1180
c          Z   =  0.5d0*Z   + 0.5*K.[D  - D  ]
 
1181
c          
 
1182
c          cp          cv           cp   pv
 
1183
c          Z   =  0.5d0*Z   - 0.5*K.[D  - D  ]
 
1184
c          
 
1185
c         Add the Fock matrices together overwriting the density
 
1186
c         matrices to form the results above
 
1187
c          
 
1188
c         Closed-Virtual bit   
 
1189
          if (debug) then
 
1190
           if (ga_nodeid().eq.0)
 
1191
     &     write(*,*) '--------- g_fck -------- START'
 
1192
           call ga_print(g_fock(ind)) 
 
1193
           if (ga_nodeid().eq.0)
 
1194
     &     write(*,*) '--------- g_fck -------- END'
 
1195
           if (ga_nodeid().eq.0)
 
1196
     &     write(*,*) '--------- g_vecs -------- START'
 
1197
           call ga_print(g_movec) 
 
1198
           if (ga_nodeid().eq.0)
 
1199
     &     write(*,*) '--------- g_vecs -------- END'
 
1200
          endif ! end-if-debug
 
1201
 
 
1202
          call ga_zero(g_tmp) ! added-04-23-12    
 
1203
          call nga_matmul_patch('n','n',four,zero,
 
1204
     &                          g_fock(ind),alo,ahi,
 
1205
     &                          g_movec    ,blo,bhi,
 
1206
     &                          g_tmp      ,blo,bhi)
 
1207
 
 
1208
          if (debug) then
 
1209
           if (ga_nodeid().eq.0)
 
1210
     &     write(*,*) '--------- FnnCno -------- START'
 
1211
           call ga_print(g_tmp) 
 
1212
           if (ga_nodeid().eq.0)
 
1213
     &     write(*,*) '--------- FnnCno -------- END'
 
1214
          endif ! end-if-debug
 
1215
 
 
1216
          call ga_zero(g_tmp1)
 
1217
          call ga_matmul_patch('t','n',one,zero,
 
1218
     $                         g_movec,voff,nmo ,1,nbf,     ! MO coefficients
 
1219
     $                         g_tmp  ,1   ,nbf ,1,nclosed, ! result from step 1 
 
1220
     $                         g_tmp1 ,1   ,nvir,1,nclosed) ! vir-occ Fock matrix
 
1221
 
 
1222
         if (debug) then
 
1223
          if (ga_nodeid().eq.0) then
 
1224
           write(*,3701) cnt,ivec,ipm
 
1225
3701     format('----- CvnFnnCno(',i3,',',i3,',',i3,')------START')
 
1226
          endif
 
1227
           call ga_print(g_tmp1)
 
1228
          if (ga_nodeid().eq.0) then
 
1229
           write(*,3702) cnt,ivec,ipm
 
1230
3702     format('----- CvnFnnCno(',i3,',',i3,',',i3,')------END')
 
1231
          endif
 
1232
         endif ! end-if-debug
 
1233
 
 
1234
          if      (cnt.eq.1) then
 
1235
 
 
1236
           if (debug) then
 
1237
            if (ga_nodeid().eq.0)        
 
1238
     &      write(*,*) '------- g_Az1-re-BEF(',ipm,')------ START' 
 
1239
            call ga_print(g_Az1)
 
1240
            if (ga_nodeid().eq.0)
 
1241
     &      write(*,*) '------- g_Az1-re-BEF(',ipm,')------ END'
 
1242
           endif ! end-if-debug
 
1243
 
 
1244
c Note.- The operation below does:
 
1245
c        g_ax_re= g_ax_re + 4 [4 C^T F C]  --> I am not sure if this is right.
 
1246
c          call ga_mat_to_vec(g_tmp1,1,nvir,1,nclosed,
 
1247
c     $                       g_ax_re(ipm),1,ivec,four,'+')
 
1248
 
 
1249
          call update_gz_reorim1(g_Az1,  ! out: = complx(g_xre,g_xim)
 
1250
     &                           g_tmp1, ! in : real      arr 
 
1251
     &                           1,      ! in  : =1 -> re =2 -> im
 
1252
     &                           nsub,   ! in : index to sub-block in g_z
 
1253
     &                           ipm,    ! in : = 1 or 2 index for component
 
1254
     &                           vlen,   ! in : = nocc*nvir
 
1255
     &                           four,   ! in  : scaling factor
 
1256
     &                           nvir,
 
1257
     &                           nclosed,
 
1258
     &                           ivec)
 
1259
 
 
1260
          if (debug) then
 
1261
           if (ga_nodeid().eq.0)        
 
1262
     &     write(*,*) '------- g_Az1-re-AFT(',ipm,')------ START' 
 
1263
           call ga_print(g_Az1)
 
1264
           if (ga_nodeid().eq.0)
 
1265
     &     write(*,*) '------- g_Az1-re-AFT(',ipm,')------ END'
 
1266
          endif ! end-if-debug
 
1267
 
 
1268
          else if (cnt.eq.2) then
 
1269
 
 
1270
          if (debug) then
 
1271
           if (ga_nodeid().eq.0)        
 
1272
     &     write(*,*) '------- g_Az1-im-BEF(',ipm,')------ START' 
 
1273
           call ga_print(g_Az1)
 
1274
           if (ga_nodeid().eq.0)
 
1275
     &     write(*,*) '------- g_Az1-im-BEF(',ipm,')------ END'
 
1276
          endif ! end-if-debug
 
1277
 
 
1278
c          call ga_mat_to_vec(g_tmp1,1,nvir,1,nclosed,
 
1279
c     $                       g_ax_im(ipm),1,ivec,four,'+')
 
1280
 
 
1281
          call update_gz_reorim1(g_Az1,  ! out: = complx(g_xre,g_xim)
 
1282
     &                           g_tmp1, ! in : real      arr 
 
1283
     &                           2,      ! in  : =1 -> re =2 -> im
 
1284
     &                           nsub,   ! in : index to sub-block in g_z
 
1285
     &                           ipm,    ! in : = 1 or 2 index for component
 
1286
     &                           vlen,   ! in : = nocc*nvir
 
1287
     &                           four,   ! in  : scaling factor
 
1288
     &                           nvir,
 
1289
     &                           nclosed,
 
1290
     &                           ivec)
 
1291
 
 
1292
          if (debug) then
 
1293
           if (ga_nodeid().eq.0)        
 
1294
     &     write(*,*) '------- g_Az1-im-AFT(',ipm,')------ START' 
 
1295
           call ga_print(g_Az1)
 
1296
           if (ga_nodeid().eq.0)
 
1297
     &     write(*,*) '------- g_Az1-im-AFT(',ipm,')------ END'
 
1298
          endif ! end-if-debug
 
1299
 
 
1300
          endif ! end-if-cnt   
 
1301
        enddo ! end-loop-ipm lopp in +/- components
 
1302
       enddo ! end-loop-ivec-loop in field directions
 
1303
      enddo ! end-loop-cnt
 
1304
      if (debug) then
 
1305
       do ipm=1,ncomp
 
1306
           if (ga_nodeid().eq.0)        
 
1307
     &     write(*,*) '------- g_Az1-1(',ipm,')------ START' 
 
1308
           call ga_print(g_Az1)
 
1309
           if (ga_nodeid().eq.0)
 
1310
     &     write(*,*) '------- g_Az1-1(',ipm,')------ END'
 
1311
       enddo ! end-loop-ipm
 
1312
      endif ! end-if-debug
 
1313
 
 
1314
c      if (debug) then
 
1315
c       if (ga_nodeid().eq.0)
 
1316
c     &  write(*,*) 'STOP-test'
 
1317
c       stop
 
1318
c      endif ! end-if-debug
 
1319
 
 
1320
      do ipm = 1,nblock
 
1321
        if (.not. ga_destroy(g_dens(ipm))) call errquit(
 
1322
     &      'rohf_hessv3: ga_destroy failed g_dens',0,GA_ERR)
 
1323
        if (.not. ga_destroy(g_fock(ipm))) call errquit
 
1324
     &     ('rohf_hessv3: ga_destroy failed g_fock',0,GA_ERR)
 
1325
      enddo ! end-loop-ipm    
 
1326
        if (.not.ga_destroy(g_tmp)) call errquit(
 
1327
     &      'rohf_hessv3: ga_destroy failed g_tmp',0,GA_ERR)
 
1328
        if (.not.ga_destroy(g_tmp1)) call errquit(
 
1329
     &      'rohf_hessv3: ga_destroy failed g_tmp',0,GA_ERR)
 
1330
      return     
 
1331
      end
 
1332
 
 
1333
      subroutine rohf_hessv_2e3_opt(
 
1334
     &                  basis,   ! in: basis handle
 
1335
     &                  geom,    ! in: geom  handle
 
1336
     &                  nbf,     ! in: nr. basis functions
 
1337
     &                  nmo,     ! in: nr. MOs
 
1338
     &                  nclosed, ! in: nr. double occupied MOs
 
1339
     &                  nopen,   ! in: nr. single occupied MOs
 
1340
     $                  g_movec, ! in: MO coefficients
 
1341
     &                  oskel,   ! in: =.true. ->
 
1342
     &                  noskew,  ! in: =.true. -> symmetric density matrix
 
1343
     &                  g_x_re,  ! in: 
 
1344
     &                  g_x_im,  ! in: 
 
1345
     &                  g_ax_re, ! in/out: Hessian product
 
1346
     &                  g_ax_im, ! in/out: Hessian product
 
1347
     &                  acc,     ! in: accuracy Fock construction
 
1348
     &                  limag,   ! in: =.true. -> imaginary component allowed
 
1349
     &                  lifetime)! in: =.true. -> RE-IM =.false -> RE
 
1350
c
 
1351
c   Purpose: Optimization of rohf_hessv_2e3()
 
1352
c   Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
1353
c   Date   : 03-28-12
 
1354
c   Note.- Modifying rohf_hessv_2e3, reducing computation of
 
1355
c          two-electron integrals by putting together 
 
1356
c          symmetric+antisymmetric perturbed densities: g_dens(ipm) ipm=1,2
 
1357
c          and doing a single call to shell_fock_build() when using
 
1358
c          the routine shell_fock_build2().
 
1359
c --> Experimental (not published yet)
 
1360
 
 
1361
      implicit none
 
1362
#include "errquit.fh"
 
1363
#include "mafdecls.fh"
 
1364
#include "global.fh"
 
1365
#include "util.fh"
 
1366
#include "cscfps.fh"
 
1367
#include "rtdb.fh"
 
1368
#include "bgj.fh"
 
1369
#include "stdio.fh"
 
1370
#include "case.fh"
 
1371
c     
 
1372
c     Return the ROHF orbital 2e-Hessian vector product, g_ax = A * g_x
 
1373
c
 
1374
c ... jochen: modified version of rohf_hessv_2e2 which keeps track
 
1375
c     of two sets of input vectors that couple via the density matrix.
 
1376
c     one could likely save some memory here by re-using temp arrays
 
1377
c ... jochen: Also made modifications to calculate imaginary terms due
 
1378
c     to finite lifetime damping   
 
1379
ccccccccccccccc This code does NOT work for open shell!!!!!ccccccccccccccccc
 
1380
      integer basis, geom       ! basis & geom handle
 
1381
      integer nclosed,nvir,nopen! Basis size and occupation
 
1382
      integer nbf,nmo           ! No. of linearly dependent MOs
 
1383
      integer g_movec           ! MO coefficients
 
1384
      logical oskel
 
1385
      integer g_x_re(2),        !
 
1386
     &        g_x_im(2)         ! Argument
 
1387
      integer g_ax_re(2),       ! Hessian product
 
1388
     &        g_ax_im(2)        ! Hessian product
 
1389
      double precision acc      ! Accuracy of "Fock" construction
 
1390
      logical limag             ! imaginary perturbation?    
 
1391
      integer voff,xoff
 
1392
      integer ivec,nvec,nmul,gtype,vlen
 
1393
      integer g_tmp,g_tmp1,g_dens(4),g_fock(4)
 
1394
      double precision tol2e
 
1395
      logical odebug,oprint,noskew
 
1396
      integer dims(3),chunk(3), 
 
1397
     &        alo(3),ahi(3), 
 
1398
     &        blo(2),bhi(2)
 
1399
      integer ga_create_atom_blocked
 
1400
      external ga_create_atom_blocked,
 
1401
     &         get_undosymm_fock,update_ax_fock,
 
1402
     &         shell_fock_build2    
 
1403
      double precision one,zero,mone,four,half,mhalf,two,mtwo
 
1404
      double precision itol_floor, itol_ceil
 
1405
      parameter(itol_floor=1.d-15, itol_ceil=1.d-3)
 
1406
      parameter (one=1.0d0, mone=-1.0d0, zero=0.0d0, four=4.0d0)
 
1407
      parameter (half=0.5d0, mhalf=-0.5d0, two=2.0d0, mtwo=-2.0d0)
 
1408
      integer ipm ! counter for density matrix components
 
1409
      character*(255) cstemp
 
1410
      integer g_pmats(2),g_pmata(2),g_h1mat(2),
 
1411
     &        cnt,ind,indx(2,2),
 
1412
     &        npol,nset,nblock,ncomp
 
1413
      double precision tenm6,coef(2,2)
 
1414
      parameter (tenm6 = 1d-6)
 
1415
      logical lifetime,debug
 
1416
      data npol /1/ ! for restricted calculations
 
1417
      data indx /1,2, ! indx(1,1),indx(1,2)
 
1418
     &           3,4/ ! indx(2,1),indx(2,2)
 
1419
      call ga_inquire(g_x_re(1),gtype,vlen,nvec) ! out: nvec,vlen
 
1420
 
 
1421
      ncomp=2 ! using two components
 
1422
      nmul=1
 
1423
      if (npol.eq. 2) nmul=2
 
1424
      nset  =1 ! for RE          
 
1425
      nblock=2 ! for RE
 
1426
      if (lifetime) then
 
1427
      nset  =2 ! for RE-IM
 
1428
      nblock=4 ! for RE-IM
 
1429
      endif
 
1430
 
 
1431
      if (oskel) 
 
1432
     $   call errquit('rohf_h2e3: no way',0, UNKNOWN_ERR)
 
1433
 
 
1434
      debug = (.false. .and. ga_nodeid().eq.0) ! for code development
 
1435
 
 
1436
c      debug=.true. ! for debugging
 
1437
 
 
1438
      if (debug) then
 
1439
       if (ga_nodeid().eq.0) then
 
1440
        write(*,2001) limag,lifetime,nset,nblock
 
1441
 2001   format('(limag,lifetime,nset,nblock)=(',
 
1442
     &           L1,',',L1,',',i3,',',i3,')')
 
1443
       endif
 
1444
      endif ! end-if-debug
 
1445
 
 
1446
      oprint= util_print('rohf_hessv2',print_debug)
 
1447
 
 
1448
      if (nopen.ne.0) call errquit
 
1449
     $     ('rohf_h2e3: does not work for open shells',nopen,
 
1450
     &       UNKNOWN_ERR)
 
1451
      odebug = util_print('rohf_hessv', print_debug)    
 
1452
      tol2e = min(max(acc,itol_floor),itol_ceil)
 
1453
      nvir = nmo - nclosed - nopen
 
1454
      voff = nclosed + nopen + 1
 
1455
 
 
1456
      dims(1)  = nbf
 
1457
      dims(2)  = nbf
 
1458
      chunk(1) = dims(1)
 
1459
      chunk(2) = -1
 
1460
c ------ create scratch GA arrays (g_pmats,g_mata,g_h1mat---START
 
1461
      do ipm = 1,ncomp
 
1462
        write(cstemp,'(a,i1)') 'pmats_',ipm
 
1463
        if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk,
 
1464
     &     g_pmats(ipm))) call 
 
1465
     &     errquit('rohf_h2e3: nga_create failed '//cstemp(1:7),
 
1466
     &     0,GA_ERR)
 
1467
        call ga_zero(g_pmats(ipm))
 
1468
        write(cstemp,'(a,i1)') 'pmata_',ipm
 
1469
        if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk,
 
1470
     &     g_pmata(ipm))) call 
 
1471
     &     errquit('rohf_h2e3: nga_create failed '//cstemp(1:7),
 
1472
     &     0,GA_ERR)
 
1473
        call ga_zero(g_pmata(ipm))
 
1474
        write(cstemp,'(a,i1)') 'h1mat_',ipm
 
1475
        if (.not.nga_create(MT_DBL,2,dims,cstemp(1:7),chunk,
 
1476
     &     g_h1mat(ipm))) call 
 
1477
     &     errquit('rohf_h2e3: nga_create failed '//cstemp(1:7),
 
1478
     &     0,GA_ERR)
 
1479
        call ga_zero(g_h1mat(ipm))
 
1480
      enddo ! end-loop-ipm
 
1481
c ------ create scratch GA arrays (g_pmats,g_mata,g_h1mat---END
 
1482
      dims(1)  = nvec
 
1483
      dims(2)  = nbf
 
1484
      dims(3)  = nbf
 
1485
      chunk(1) = dims(1)
 
1486
      chunk(2) = -1
 
1487
      chunk(3) = -1     
 
1488
      do ipm = 1,nblock ! =2 or 4
 
1489
c ... allocate g_dens=[g_dens_re g_dens_im]
 
1490
        if (.not. nga_create (MT_DBL,3,dims,'CPKS dens',chunk,
 
1491
     &     g_dens(ipm)))
 
1492
     &     call errquit('rohf_h2e3: could not allocate g_dens',555,
 
1493
     &     GA_ERR)
 
1494
        call ga_zero(g_dens(ipm))
 
1495
c ... allocate g_fock=[g_fock_re g_fock_im]
 
1496
        if (.not. nga_create (MT_DBL,3,dims,'Fockv',chunk,
 
1497
     &     g_fock(ipm)))
 
1498
     &     call errquit('rohf_h2e3: could not allocate g_fock',555,
 
1499
     &     GA_ERR)
 
1500
        call ga_zero(g_fock(ipm))
 
1501
      enddo ! end-loop-ipm
 
1502
 
 
1503
        if (debug) then
 
1504
         if (ga_nodeid().eq.0)
 
1505
     &    write(*,*) 'BEF get_dens_reorim-RE'
 
1506
        endif ! end-if-debug
 
1507
 
 
1508
      call get_dens_reorim(
 
1509
     &                    g_dens, ! in/ou: perturbed density matrix
 
1510
     &                    1,      ! in   : =1 1st block RE
 
1511
     &                    g_x_re, ! in   : RE
 
1512
     &                    g_movec,! in   : MO coefficients
 
1513
     &                    nbf,    ! in   : nr. basis functions
 
1514
     &                    nmo,    ! in   : nr. MOs
 
1515
     &                    1,      ! in   : for ipol=1 -> istart=1
 
1516
     &                    nclosed,! in   : nr. occupied MOs
 
1517
     &                    nvir,   ! in   : nr. virtual  MOs 
 
1518
     &                    nvec,   ! in   : nr. directions (x,y,z)
 
1519
     &                    npol,   ! in   : nr. polarizations
 
1520
     &                    limag,  ! in   : = .true. imaginary allowed   
 
1521
     &                    g_pmats,! in   : scratch GA array
 
1522
     &                    g_pmata,! in   : scratch GA array   
 
1523
     &                    g_h1mat)! in   : scratch GA array   
 
1524
 
 
1525
         if (debug) then
 
1526
          if (ga_nodeid().eq.0)
 
1527
     &     write(*,*) 'BEF get_dens_reorim-IM'
 
1528
         endif ! end-if-debug
 
1529
 
 
1530
      if (lifetime) then
 
1531
 
 
1532
       call get_dens_reorim(
 
1533
     &                    g_dens, ! in/ou: perturbed density matrix
 
1534
     &                    2,      ! in   : =2 2nd block IM
 
1535
     &                    g_x_im, ! in   : IM
 
1536
     &                    g_movec,! in   : MO coefficients
 
1537
     &                    nbf,    ! in   : nr. basis functions
 
1538
     &                    nmo,    ! in   : nr. MOs
 
1539
     &                    1,      ! in   : for ipol=1 -> istart=1
 
1540
     &                    nclosed,! in   : nr. occupied MOs
 
1541
     &                    nvir,   ! in   : nr. virtual  MOs 
 
1542
     &                    nvec,   ! in   : nr. directions (x,y,z)
 
1543
     &                    npol,   ! in   : nr. polarizations
 
1544
     &                    limag,  ! in   : = .true. imaginary allowed   
 
1545
     &                    g_pmats,! in   : scratch GA array
 
1546
     &                    g_pmata,! in   : scratch GA array   
 
1547
     &                    g_h1mat)! in   : scratch GA array   
 
1548
 
 
1549
      endif ! end-if-lifetime
 
1550
 
 
1551
c         if (ga_nodeid().eq.0) 
 
1552
c     &     write(*,*) 'FA-BEF shell_fock_build2 ...'
 
1553
 
 
1554
      if (debug) then
 
1555
        do ipm=1,nblock
 
1556
         if (ga_nodeid().eq.0)
 
1557
     &    write(*,*) '--------g_dens(',ipm,')-------START'
 
1558
         call ga_print(g_dens(ipm))
 
1559
         if (ga_nodeid().eq.0)
 
1560
     &    write(*,*) '--------g_dens(',ipm,')-------END'
 
1561
        enddo ! end-loop-ipm
 
1562
      endif ! end-if-debug
 
1563
 
 
1564
      call shell_fock_build2(g_fock,  ! out: Fock    matrices
 
1565
     &                       g_dens,  ! in : density matrices
 
1566
     &                       geom,    ! in : geom  handle
 
1567
     &                       basis,   ! in : basis handle
 
1568
     &                       nbf,     ! in : nr. basis functions
 
1569
     &                       nvec,    ! in : nr. vecs (x,y,z)
 
1570
     &                       npol,    ! in : npol=1 for R-DFT =2 for U-DFT
 
1571
     &                       ncomp,   ! in : nr. components
 
1572
     &                       nblock,  ! in : nr. of g_dens,g_fock blocks
 
1573
     &                       .true.,  ! in : =.true. for symm dens
 
1574
     &                       tol2e,   ! in :
 
1575
     &                       debug)   ! in : = .true. -> debugging printouts
 
1576
 
 
1577
c       if (ga_nodeid().eq.0) 
 
1578
c     &  write(*,*) 'FA-AFT shell_fock_build2 ...'
 
1579
 
 
1580
      if (debug) then
 
1581
       do ipm=1,nblock
 
1582
         if (ga_nodeid().eq.0)        
 
1583
     &    write(*,*) '------- g_fock-0(',ipm,')------ START' 
 
1584
         call ga_print(g_fock(ipm))
 
1585
         if (ga_nodeid().eq.0)
 
1586
     &    write(*,*) '------- g_fock-0(',ipm,')------ END'
 
1587
       enddo ! end-loop-ipm
 
1588
      endif ! end-if-debug
 
1589
 
 
1590
      call get_undosymm_fock(
 
1591
     &            g_fock,  ! in/ou: fock matrix
 
1592
     &            nset,    ! in   : =1 g_x is real, =2 g_x is complex (g_x_re,g_x_im)
 
1593
     &            nvec,    ! in   : nr. directions (x,y,z)
 
1594
     &            nbf,     ! in   : nr. basis functions
 
1595
     &            npol,    ! in   : nr. polarizations
 
1596
     &            nmul,    ! in   : =1 npol=1 =2 npol=2 (acc. JK terms)
 
1597
     &            g_pmats, ! in   : scratch GA array
 
1598
     &            limag)   ! in   : =.true. imaginary comp. exists
 
1599
 
 
1600
c ------- Remove GA arrays:
 
1601
      do ipm = 1,ncomp
 
1602
       if (.not.ga_destroy(g_pmats(ipm))) call errquit(
 
1603
     &     'rohf_hessv3: ga_destroy failed g_pmats',0,GA_ERR)       
 
1604
       if (.not.ga_destroy(g_pmata(ipm))) call errquit(
 
1605
     &      'rohf_hessv3: ga_destroy failed g_pmata',0,GA_ERR)
 
1606
       if (.not.ga_destroy(g_h1mat(ipm)))  call errquit(
 
1607
     &      'rohf_hessv3: ga_destroy failed g_h1mat',0,GA_ERR)
 
1608
      enddo ! end-loop-ipm
 
1609
 
 
1610
      if (debug) then
 
1611
       do ipm=1,nblock
 
1612
         if (ga_nodeid().eq.0)        
 
1613
     &    write(*,*) '------- g_fock-unsym(',ipm,')------ START' 
 
1614
         call ga_print(g_fock(ipm))
 
1615
         if (ga_nodeid().eq.0)
 
1616
     &    write(*,*) '------- g_fock-unsym(',ipm,')------ END'
 
1617
       enddo ! end-loop-ipm
 
1618
      endif ! end-if-debug
 
1619
        
 
1620
      if (debug) then
 
1621
       do ipm=1,ncomp
 
1622
         if (ga_nodeid().eq.0)        
 
1623
     &    write(*,*) '------- g_ax_re-0(',ipm,')------ START' 
 
1624
         call ga_print(g_ax_re(ipm))
 
1625
         if (ga_nodeid().eq.0)
 
1626
     &    write(*,*) '------- g_ax_re-0(',ipm,')------ END'
 
1627
       enddo ! end-loop-ipm
 
1628
       if (lifetime) then
 
1629
       do ipm=1,ncomp
 
1630
         if (ga_nodeid().eq.0)        
 
1631
     &    write(*,*) '------- g_ax_im-0(',ipm,')------ START' 
 
1632
         call ga_print(g_ax_im(ipm))
 
1633
         if (ga_nodeid().eq.0)
 
1634
     &    write(*,*) '------- g_ax_im-0(',ipm,')------ END'
 
1635
       enddo ! end-loop-ipm
 
1636
       endif ! end-if-lifetime
 
1637
      endif ! end-if-debug
 
1638
c     
 
1639
c     start loop over components of perturbing field
 
1640
c   
 
1641
      g_tmp = ga_create_atom_blocked(geom, basis,'rohf_h2e3: tmp')
 
1642
      g_tmp1= ga_create_atom_blocked(geom, basis,'rohf_h2e3: tmp1')
 
1643
      alo(2) = 1
 
1644
      ahi(2) = nbf
 
1645
      alo(3) = 1
 
1646
      ahi(3) = nbf
 
1647
      blo(1) = 1
 
1648
      bhi(1) = nbf
 
1649
      blo(2) = 1
 
1650
      bhi(2) = nclosed  
 
1651
      do cnt=1,nset
 
1652
       do ivec = 1, nvec
 
1653
        alo(1) = ivec
 
1654
        ahi(1) = ivec      
 
1655
        do ipm = 1,ncomp ! loop over Fock matrix components +/- here     
 
1656
          ind=indx(ipm,cnt)
 
1657
 
 
1658
           if (debug) then
 
1659
            if (ga_nodeid().eq.0) then
 
1660
             write(*,117) cnt,ivec,ipm,ind
 
1661
 117         format('XX:(cnt,ivec,ipm,ind)=(',
 
1662
     &              i3,',',i3,',',i3,',',i3,')')
 
1663
            endif
 
1664
           endif ! end-if-debug
 
1665
c         
 
1666
c          P      =  4(ij|kl) - (ik|jl) - (il|kj)
 
1667
c           ij,kl
 
1668
c     
 
1669
c          K      =  (ik|jl) + (il|kj)
 
1670
c           ij,kl
 
1671
c     
 
1672
c          cv         cv          pv   cp
 
1673
c          Z   =  2P.[D  ]  +  P.[D  + D  ]
 
1674
c     
 
1675
c          pv          cv           cp   pv
 
1676
c          Z   =  0.5d0*Z   + 0.5*K.[D  - D  ]
 
1677
c          
 
1678
c          cp          cv           cp   pv
 
1679
c          Z   =  0.5d0*Z   - 0.5*K.[D  - D  ]
 
1680
c          
 
1681
c         Add the Fock matrices together overwriting the density
 
1682
c         matrices to form the results above
 
1683
c          
 
1684
c         Closed-Virtual bit   
 
1685
        if (debug) then
 
1686
         if (ga_nodeid().eq.0)
 
1687
     &    write(*,*) '--------- g_fck -------- START'
 
1688
         call ga_print(g_fock(ind)) 
 
1689
         if (ga_nodeid().eq.0)
 
1690
     &    write(*,*) '--------- g_fck -------- END'
 
1691
         if (ga_nodeid().eq.0)
 
1692
     &    write(*,*) '--------- g_vecs -------- START'
 
1693
         call ga_print(g_movec) 
 
1694
         if (ga_nodeid().eq.0)
 
1695
     &    write(*,*) '--------- g_vecs -------- END'
 
1696
        endif ! end-if-debug
 
1697
          call ga_zero(g_tmp) ! added-04-23-12         
 
1698
          call nga_matmul_patch('n','n',four,zero,
 
1699
     &                          g_fock(ind),alo,ahi,
 
1700
     &                          g_movec    ,blo,bhi,
 
1701
     &                          g_tmp      ,blo,bhi)
 
1702
        if (debug) then
 
1703
         if (ga_nodeid().eq.0)
 
1704
     &    write(*,*) '--------- FnnCno -------- START'
 
1705
         call ga_print(g_tmp) 
 
1706
         if (ga_nodeid().eq.0)
 
1707
     &    write(*,*) '--------- FnnCno -------- END'
 
1708
        endif ! end-if-debug
 
1709
 
 
1710
          call ga_zero(g_tmp1)
 
1711
          call ga_matmul_patch('t','n',one,zero,
 
1712
     $                         g_movec,voff,nmo ,1,nbf,     ! MO coefficients
 
1713
     $                         g_tmp  ,1   ,nbf ,1,nclosed, ! result from step 1 
 
1714
     $                         g_tmp1 ,1   ,nvir,1,nclosed) ! vir-occ Fock matrix
 
1715
 
 
1716
         if (debug) then
 
1717
          if (ga_nodeid().eq.0) then
 
1718
           write(*,3701) cnt,ivec,ipm
 
1719
3701     format('----- CvnFnnCno(',i3,',',i3,',',i3,')------START')
 
1720
          endif
 
1721
           call ga_print(g_tmp1)
 
1722
          if (ga_nodeid().eq.0) then
 
1723
           write(*,3702) cnt,ivec,ipm
 
1724
3702     format('----- CvnFnnCno(',i3,',',i3,',',i3,')------END')
 
1725
          endif
 
1726
         endif ! end-if-debug
 
1727
 
 
1728
          if      (cnt.eq.1) then
 
1729
 
 
1730
          if (debug) then
 
1731
           if (ga_nodeid().eq.0)
 
1732
     &     write(*,*) '--------- g_ax-re-BEF-------- START'
 
1733
           call ga_print(g_ax_re(ipm)) 
 
1734
           if (ga_nodeid().eq.0)
 
1735
     &     write(*,*) '--------- g_ax-re-BEF-------- END'
 
1736
          endif ! end-if-debug
 
1737
 
 
1738
c Note.- The operation below does:
 
1739
c        g_ax_re= g_ax_re + 4 [4 C^T F C]  --> I am not sure if this is right.
 
1740
          call ga_mat_to_vec(g_tmp1,1,nvir,1,nclosed,
 
1741
     $                       g_ax_re(ipm),1,ivec,four,'+')
 
1742
 
 
1743
          if (debug) then
 
1744
           if (ga_nodeid().eq.0)
 
1745
     &      write(*,*) '--------- g_ax-re-AFT-------- START'
 
1746
           call ga_print(g_ax_re(ipm)) 
 
1747
           if (ga_nodeid().eq.0)
 
1748
     &      write(*,*) '--------- g_ax-re-AFT-------- END'
 
1749
          endif ! end-if-debug
 
1750
 
 
1751
          else if (cnt.eq.2) then
 
1752
 
 
1753
          if (debug) then
 
1754
           if (ga_nodeid().eq.0)
 
1755
     &      write(*,*) '--------- g_ax-im-BEF-------- START'
 
1756
           call ga_print(g_ax_im(ipm)) 
 
1757
           if (ga_nodeid().eq.0)
 
1758
     &      write(*,*) '--------- g_ax-im-BEF-------- END'
 
1759
          endif ! end-if-debug
 
1760
 
 
1761
          call ga_mat_to_vec(g_tmp1,1,nvir,1,nclosed,
 
1762
     $                       g_ax_im(ipm),1,ivec,four,'+')
 
1763
 
 
1764
          if (debug) then
 
1765
           if (ga_nodeid().eq.0)
 
1766
     &      write(*,*) '--------- g_ax-im-AFT-------- START'
 
1767
           call ga_print(g_ax_im(ipm)) 
 
1768
           if (ga_nodeid().eq.0)
 
1769
     &      write(*,*) '--------- g_ax-im-AFT-------- END'
 
1770
          endif ! end-if-debug
 
1771
 
 
1772
          endif ! end-if-cnt   
 
1773
        enddo ! end-loop-ipm lopp in +/- components
 
1774
       enddo ! end-loop-ivec-loop in field directions
 
1775
      enddo ! end-loop-cnt
 
1776
 
 
1777
      if (debug) then
 
1778
       do ipm=1,ncomp
 
1779
         if (ga_nodeid().eq.0)        
 
1780
     &    write(*,*) '------- g_ax_re-1(',ipm,')------ START' 
 
1781
         call ga_print(g_ax_re(ipm))
 
1782
         if (ga_nodeid().eq.0)
 
1783
     &    write(*,*) '------- g_ax_re-1(',ipm,')------ END'
 
1784
       enddo ! end-loop-ipm
 
1785
       if (lifetime) then
 
1786
       do ipm=1,ncomp
 
1787
         if (ga_nodeid().eq.0)        
 
1788
     &    write(*,*) '------- g_ax_im-1(',ipm,')------ START' 
 
1789
         call ga_print(g_ax_im(ipm))
 
1790
         if (ga_nodeid().eq.0)
 
1791
     &    write(*,*) '------- g_ax_im-1(',ipm,')------ END'
 
1792
       enddo ! end-loop-ipm
 
1793
       endif ! end-if-lifetime
 
1794
      endif ! end-if-debug
 
1795
 
 
1796
c      if (debug) then
 
1797
c       if (ga_nodeid().eq.0)
 
1798
c     &  write(*,*) 'STOP-test'
 
1799
c       stop
 
1800
c      endif ! end-if-debug
 
1801
 
 
1802
      do ipm = 1,nblock
 
1803
        if (.not. ga_destroy(g_dens(ipm))) call errquit(
 
1804
     &      'rohf_hessv3: ga_destroy failed g_dens',0,GA_ERR)
 
1805
        if (.not. ga_destroy(g_fock(ipm))) call errquit
 
1806
     &     ('rohf_hessv3: ga_destroy failed g_fock',0,GA_ERR)
 
1807
      enddo ! end-loop-ipm    
 
1808
        if (.not.ga_destroy(g_tmp)) call errquit(
 
1809
     &      'rohf_hessv3: ga_destroy failed g_tmp',0,GA_ERR)
 
1810
        if (.not.ga_destroy(g_tmp1)) call errquit(
 
1811
     &      'rohf_hessv3: ga_destroy failed g_tmp',0,GA_ERR)
 
1812
      return     
 
1813
      end
 
1814
 
 
1815
      subroutine get_dens_reorim(
 
1816
     &                    g_dens, ! out  : perturbed density matrix
 
1817
     &                    cnt,    ! in/ou: counter of g_dens, =1 or 2
 
1818
     &                    g_x,    ! in   : 
 
1819
     &                    g_movec,! in   : MO coefficients
 
1820
     &                    nbf,    ! in   : nr. basis functions
 
1821
     &                    nmo,    ! in   : nr. MOs
 
1822
     &                    istart, ! in   : shift nocc-nvirt block
 
1823
     &                    nocc,   ! in   : nr. occupied MOs
 
1824
     &                    nvir,   ! in   : nr. virtual  MOs 
 
1825
     &                    nvec,   ! in   : nr. directions (x,y,z)
 
1826
     &                    ipol,   ! in   : nr. polarizations
 
1827
     &                    limag,  ! in   : = .true. imaginary allowed   
 
1828
     &                    g_pmats,! in   : scratch GA array
 
1829
     &                    g_pmata,! in   : scratch GA array - NOT USED
 
1830
     &                    g_h1mat)! in   : scratch GA array   
 
1831
c
 
1832
c Purpose: Calculate perturbed density matrix
 
1833
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
1834
c Date   : 03-15-12
 
1835
c --> Experimental (not published yet)
 
1836
 
 
1837
      implicit none
 
1838
#include "errquit.fh"
 
1839
#include "mafdecls.fh"
 
1840
#include "global.fh"
 
1841
#include "util.fh"
 
1842
#include "cscfps.fh"
 
1843
#include "rtdb.fh"
 
1844
#include "bgj.fh"
 
1845
#include "stdio.fh"
 
1846
#include "case.fh"
 
1847
      integer g_x(2)            ! Argument: g_x_re or g_x_im
 
1848
      integer g_dens(*)         ! size= 2 RE only or 4 RE+IM
 
1849
      integer nbf,nmo,nocc,nvir
 
1850
      integer g_movec           ! MO coefficients
 
1851
      integer dims(3),chunk(3), 
 
1852
     &        alo(3),ahi(3), 
 
1853
     &        blo(2),bhi(2)
 
1854
      integer ivec,nvec,ipm,ncomp,
 
1855
     &        ipol, ! =1 for Alpha =2 for Beta
 
1856
     &        shift,cnt,ind,indx(2,2)
 
1857
      character*(255) cstemp
 
1858
      integer g_pmats(2),g_pmata(2),g_h1mat(2)
 
1859
      integer istart,iend
 
1860
      double precision tenm6,coef(2,2)
 
1861
      parameter (tenm6 = 1d-6)
 
1862
      logical limag,debug
 
1863
      double precision one, zero, mone, 
 
1864
     &                 four, half, mhalf, two, mtwo
 
1865
      data indx /1,2, ! indx(1,1),indx(2,1)
 
1866
     &           3,4/ ! indx(1,2),indx(2,2) 
 
1867
      data ncomp/2/
 
1868
      parameter (one=1.0d0, mone=-1.0d0, zero=0.0d0, four=4.0d0)
 
1869
      parameter (half=0.5d0, mhalf=-0.5d0, two=2.0d0, mtwo=-2.0d0)
 
1870
      external ga_vec_to_mat,
 
1871
     &         CalcPerturbedTDPmat1_opt
 
1872
 
 
1873
c      debug=.true. ! allow debugging printouts
 
1874
      debug=.false. ! allow debugging printouts
 
1875
 
 
1876
c ----- Construct coeffs for P(S),P(A) ------- START
 
1877
      coef(1,1)= 0.5d0
 
1878
      coef(1,2)= 0.5d0
 
1879
      coef(2,1)=-0.5d0
 
1880
      coef(2,2)= 0.5d0
 
1881
      if (limag) then
 
1882
      coef(1,1)= 0.5d0
 
1883
      coef(1,2)=-0.5d0
 
1884
      coef(2,1)=-0.5d0
 
1885
      coef(2,2)=-0.5d0
 
1886
      endif ! end-if-limag
 
1887
      if (debug) then
 
1888
       if (ga_nodeid().eq.0) then
 
1889
        write(*,10) limag,
 
1890
     &              coef(1,1),coef(1,2),
 
1891
     &              coef(2,1),coef(2,2)
 
1892
 10     format('(limag,coeff)=(',L1,',',f15.8,',',
 
1893
     &         f15.8,',',f15.8,',',f15.8,')')
 
1894
       endif 
 
1895
      endif
 
1896
c ----- Construct coeffs for P(S),P(A) ------- END
 
1897
      alo(2) = 1
 
1898
      ahi(2) = nbf
 
1899
      alo(3) = 1
 
1900
      ahi(3) = nbf
 
1901
      blo(1) = 1
 
1902
      bhi(1) = nbf
 
1903
      blo(2) = 1
 
1904
      bhi(2) = nbf
 
1905
      shift=(ipol-1)*nvec
 
1906
      iend = istart + nocc*nvir - 1     
 
1907
      do ivec = 1, nvec
 
1908
c       
 
1909
c       Compute CV, PV & CP "densities" from argument vector
 
1910
c       
 
1911
c ... jochen: skip this part and place a subroutine call instead.
 
1912
c       it calculates the perturbed density matrix in the AO basis.
 
1913
c       I keep this source code here for reference; it is left
 
1914
c       unmodified from the version of rohf_hessv2 that this
 
1915
c       subroutine was created from. 
 
1916
        do ipm = 1,ncomp
 
1917
          call ga_zero(g_h1mat(ipm))
 
1918
          call ga_copy_patch('n', ! Reshape vector into matrix 
 
1919
     $                       g_x(ipm)    ,istart,iend,ivec,ivec,
 
1920
     $                       g_h1mat(ipm),1     ,nvir,1   ,nocc)
 
1921
 
 
1922
        enddo ! end-loop-ipm
 
1923
 
 
1924
       if (debug) then
 
1925
       do ipm=1,2
 
1926
        if (ga_nodeid().eq.0)
 
1927
     &   write(*,*) '----------g_h1mat(',ipm,')-----START'
 
1928
        call ga_print(g_h1mat(ipm))
 
1929
        if (ga_nodeid().eq.0)
 
1930
     &   write(*,*) '----------g_h1mat(',ipm,')-----END'
 
1931
       enddo ! end-loop-ip
 
1932
       endif ! end-if-debug
 
1933
 
 
1934
        if (debug) then
 
1935
         if (ga_nodeid().eq.0)
 
1936
     &   write(*,*) 'In rohf_hessv_2e3: BEF CalcPerturbedTDPmat1'
 
1937
         if (ga_nodeid().eq.0) then
 
1938
          write(*,667) 2,nbf,nocc,nvir,nmo,limag
 
1939
 667      format('(ncomp,nbf,nclosed,nvir,nmo,limag)=(',
 
1940
     &          i3,',',i3,',',i3,',',i3,',',i3,',',L1,')')
 
1941
         endif
 
1942
        endif ! end-if -debug
 
1943
 
 
1944
c        if (ga_nodeid().eq.0)
 
1945
c     &   write(*,*) 'Calling CalcPerturbedTDPmat1_opt!'
 
1946
 
 
1947
              call CalcPerturbedTDPmat1_opt(
 
1948
     &                 2,        ! in : nr. components to calculate
 
1949
     &                 g_pmats,  ! out: density matrix      symmetrized
 
1950
     &                 g_pmata,  ! out: density matrix  antisymmetrized
 
1951
     &                 g_h1mat,  ! in :   perturbed MO coefficients
 
1952
     &                 g_movec,  ! in : unperturbed MO coefficients
 
1953
     &                 nbf,      ! in : nr. AOs
 
1954
     &                 nocc,     ! in : nr. occupied MOs
 
1955
     &                 nvir,     ! in : nr. virtual  MOs
 
1956
     &                 nmo,      ! in : nr. MOs
 
1957
     &                 .false.,  ! in : = .true. calc. (symm,antisymm)=(pmats,pmata)
 
1958
     &                 .false.,  ! in : = .true. static response, dynamic otherwise
 
1959
     &                  limag,   ! in : = .true. if amat is imaginary instead of real
 
1960
     &                 .false.)  ! in : = .true. if amat contains occ-occ
 
1961
        
 
1962
c        if (ga_nodeid().eq.0)
 
1963
c     &   write(*,*) 'In rohf_hessv_2e3: AFT CalcPerturbedTDPmat1'
 
1964
 
 
1965
         if (debug) then
 
1966
             if (ga_nodeid().eq.0)
 
1967
     &       write(*,*) '---- g_pmats-1-------- START'
 
1968
             call ga_print(g_pmats(1))
 
1969
            if (ga_nodeid().eq.0)
 
1970
     &       write(*,*) '---- g_pmats-1-------- END'    
 
1971
c             if (ga_nodeid().eq.0)
 
1972
c     &       write(*,*) '---- g_pmata-1-------- START'
 
1973
c             call ga_print(g_pmata(1))
 
1974
c            if (ga_nodeid().eq.0)
 
1975
c     &       write(*,*) '---- g_pmata-1-------- END'   
 
1976
   
 
1977
             if (ga_nodeid().eq.0)
 
1978
     &       write(*,*) '---- g_pmats-2-------- START'
 
1979
             call ga_print(g_pmats(2))
 
1980
            if (ga_nodeid().eq.0)
 
1981
     &       write(*,*) '---- g_pmats-2-------- END'    
 
1982
c             if (ga_nodeid().eq.0)
 
1983
c     &       write(*,*) '---- g_pmata-2-------- START'
 
1984
c             call ga_print(g_pmata(2))
 
1985
c            if (ga_nodeid().eq.0)
 
1986
c     &       write(*,*) '---- g_pmata-2-------- END'        
 
1987
         endif ! end-if-debug
 
1988
 
 
1989
c next 2 lines for debugging only, to force uncoupled CPKS
 
1990
c        call ga_zero(g_pmata(1))
 
1991
c        call ga_zero(g_pmata(2))
 
1992
c       
 
1993
c       Instead of P(+) and P(-) which are both non-symmetric for
 
1994
c       non-zero frequency
 
1995
c       we will work with a symmetrized (S) and an antisymmetrized (A)
 
1996
c       component, calculate F(S) and F(A), respectively, and construct
 
1997
c       the Fock operators F(+/-) afterwards from F(S) +/- F(A).
 
1998
c       If it works for the skew-symmetric density matrix of NMR then
 
1999
c       it should work for this problem here, too
 
2000
c       note: here is one of those ominous scalings by 1/4
 
2001
c       needed to get the correct final results
 
2002
        call ga_scale(g_pmats(1),0.25d0)
 
2003
        call ga_scale(g_pmats(2),0.25d0)     
 
2004
        alo(1) = shift+ivec
 
2005
        ahi(1) = alo(1)
 
2006
c       we need to take care here of the symmetry of the density
 
2007
c       matrices depending on whether the perturbation is real
 
2008
c       or purely imaginary.
 
2009
c
 
2010
c       this works for real, symmetric, perturbations
 
2011
c       calculate P(S) = [ P(+) + P(-)]/2
 
2012
c       calculate P(A) = [-P(+) + P(-)]/2  (wrong results
 
2013
c                                          with opposite sign ...)
 
2014
         do ipm=1,ncomp
 
2015
          ind=indx(ipm,cnt)
 
2016
 
 
2017
          if (debug) then
 
2018
          if (ga_nodeid().eq.0) then
 
2019
           write(*,11) cnt,ipm,ind,ivec,coef(ipm,1),coef(ipm,2)
 
2020
 11        format('check-ind: (cnt,ipm,ind,ivec)=(',
 
2021
     &            i3,',',i3,',',i3,',',i3,')',
 
2022
     &            'coeff12=(',f15.8,',',f15.8,')')
 
2023
          endif
 
2024
          endif ! end-if-debug
 
2025
 
 
2026
          call nga_add_patch(coef(ipm,1),g_pmats(1) ,blo,bhi,
 
2027
     &                       coef(ipm,2),g_pmats(2) ,blo,bhi,
 
2028
     &                                   g_dens(ind),alo,ahi)
 
2029
          if (debug) then
 
2030
            if (ga_nodeid().eq.0) then
 
2031
             write(*,*) '---g_dens-acc(',ind,',',ivec,')-------START'
 
2032
            endif
 
2033
            call ga_print(g_dens(ind))
 
2034
            if (ga_nodeid().eq.0)
 
2035
     &      write(*,*) '----g_dens-acc(',ind,',',ivec,')-------END'
 
2036
          endif ! end-if-debug
 
2037
 
 
2038
         enddo ! end-loop-ipm   
 
2039
      enddo ! end-loop-ivec 
 
2040
      return
 
2041
      end