← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/property/giao_b1_movecs.F

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

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/property/giao_b1_movecs.F
1
 
      subroutine giao_b1_movecs(rtdb,basis,geom, g_vecB1, g_sket1)
2
 
 
3
 
c $Id: giao_b1_movecs.F 21176 2011-10-10 06:35:49Z d3y133 $
 
1
      subroutine giao_b1_movecs(rtdb,basis,geom,  ! IN
 
2
     &                          g_vecB1, g_sket1) ! OUT
 
3
c $Id: giao_b1_movecs.F 23263 2012-12-09 18:38:17Z niri $
4
4
 
5
5
c     This routine is a modification of hnd_giaox.F. Instead of
6
6
c     computing NMR shielding tensors we simply calculate the MO vectors
12
12
 
13
13
c     Note: integrals have to be initialized by the calling routine.
14
14
c     Note: the CPHF call below terminates the integrals
 
15
c
 
16
c  Written by J. Autschbach, SUNY Buffalo
 
17
c  Extension to spin-unrestricted case 
 
18
c          by F. Aquino,     Northwestern University 
 
19
c          03-15-12
 
20
c --> Experimental (not published yet)
15
21
 
16
22
      implicit none
17
 
 
18
23
#include "errquit.fh"
19
24
#include "global.fh"
20
25
#include "mafdecls.fh"
27
32
#include "prop.fh"
28
33
#include "bgj.fh"
29
34
#include "case.fh"
30
 
c
31
 
      integer rtdb    ! [in] rtdb handle
32
 
      integer basis   ! [in] basis handle
33
 
      integer geom    ! [in] geometry handle
34
 
      integer g_vecB1 ! [out] B-field perturbed MO coefficients GIAO
 
35
      integer rtdb       ! [in] rtdb handle
 
36
      integer basis      ! [in] basis handle
 
37
      integer geom       ! [in] geometry handle
 
38
      integer g_vecB1(2) ! [out] B-field perturbed MO coefficients GIAO
35
39
      integer g_sket1 ! [out] GIAO right hand side overlap derivative
36
 
 
37
40
      integer nclosed(2), nopen(2), nvirt(2), ndens, nbf, nmo
38
41
      integer ixy, ix, iy, iatom, iocc, ifld, ioff
39
42
      integer alo(3), ahi(3), blo(3), bhi(3), clo(3), chi(3)
40
43
      integer dlo(3), dhi(3)
41
44
      integer l_occ, k_occ, l_eval, k_eval
42
 
 
43
 
      integer g_dens(3), g_s10, g_d1, g_rhs, g_fock, g_u
 
45
      integer g_dens(3), g_s10, g_d1, g_rhs, g_fock, g_u(2)
44
46
      integer vectors(2), geomnew, i, j, ij, g_xc(3)
 
47
      integer vectors1(2)
 
48
      integer ndata,m1,m2
45
49
      double precision atn, tol2e, val
46
 
      double precision jfac(3),kfac(3),a(6),xfac
 
50
      double precision jfac(12),kfac(12),a(6),xfac
47
51
      character*3 scftyp
48
52
      character*16 tag
49
53
      character*32 element
50
54
      character*256 cphf_rhs, cphf_sol
51
55
      character*2 symbol
52
 
 
53
 
      double precision origin(3)
54
 
      data origin/0d0,0d0,0d0/
55
 
 
56
 
      integer nat
57
 
      parameter (nat=1)
58
 
 
 
56
c ===== added for unrestricted calc ===== START
 
57
      integer ndir,    ! nr directions (x,y,z)
 
58
     &        ntot,    ! sum_{i=1,npol} nocc(i)*nvirt(i)
 
59
     &        ispin,disp,shift,nind_jk,
 
60
     &        nocc(2), ! store nr occupations 
 
61
     &        npol     ! nr of polarizations =1 (  restricted) 
 
62
                       !                     =2 (unrestricted) calc    
 
63
      external get_d1_giao_b1,update_rhs_shfock,
 
64
     &         get_fock2e,update_rhs_threeAOints,get_vecB1,
 
65
     &         get_nocc,update_rhs_eS10,get_vecB1_opt2
 
66
c ===== added for unrestricted calc ===== END
59
67
      integer nbq, nextbq, ncosbq
60
 
 
 
68
      integer g_rhs_im,read_grhs_giaob1,n_data
 
69
      character*255 aorespfilename
 
70
      logical dft_CPHF1_read,dft_CPHF1_write
 
71
      character*(*) lbl_cphfaoresp
 
72
      parameter(lbl_cphfaoresp='aoresp_giao_b1')
61
73
      logical  cphf2, file_write_ga, file_read_ga, cphf
62
74
      external cphf2, file_write_ga, file_read_ga, cphf
63
75
 
64
 
      logical     oskel, status, debug
 
76
      logical  oskel, status, debug
65
77
      logical  xc_gotxc
66
78
      external xc_gotxc
67
79
      double precision ppm
68
80
      data tol2e   /1.0d-10/
69
 
 
70
81
c     ==================================================================
71
 
c
 
82
      ndir=3 ! nr directions (x,y,z)
72
83
      if (ga_nodeid().eq.0) write(luout,9999)
73
 
 
74
 
      call ga_zero(g_vecB1) ! initialize array for main result
75
 
 
76
84
      debug = .false. .and. (ga_nodeid().eq.0) ! special debugging
 
85
      
 
86
c      debug=.true.
77
87
 
78
88
      if (debug) then
79
89
        write (luout,*) 'giao_b1_movecs: xc_gotxc =',xc_gotxc()
80
90
        write (luout,*) 'giao_b1_movecs: use_theory =',use_theory
81
91
      end if
82
 
 
83
92
c     there is a possibility that this routine is called from tddft
84
93
c     in which case use_theory is not set. We set it to 'dft' in that case,
85
94
c     assuming that we are indeed calling from some DFT response code
89
98
     &     'giao_b1_movecs: assuming DFT/TDDFT'
90
99
        use_theory = 'dft'
91
100
      end if
92
 
 
93
101
c
94
102
c     Current CPHF does not handle symmetry 
95
103
c     Making C1 geometry and store it on rtdb
97
105
      oskel = .false.
98
106
c
99
107
c     If DFT get part of the exact exchange defined
100
 
c
101
108
      xfac = 1.0d0
102
109
      if (use_theory.eq.'dft') xfac = bgj_kfac()
 
110
      nind_jk=12
 
111
      do ifld = 1,nind_jk
 
112
        jfac(ifld) =  0.0d0       ! used in update_rhs_shfock()
 
113
        kfac(ifld) = -1.0d0*xfac  ! used in update_rhs_shfock()
 
114
c        if (ga_nodeid().eq.0) then
 
115
c         write(*,144) ifld,jfac(ifld),kfac(ifld)
 
116
c  144    format('(j,k)(',i3,')=(',f15.8,',',f15.8,')')
 
117
c        endif
 
118
      enddo
103
119
c
104
120
c     Integral initialization
105
 
c
106
 
c      call int_init(rtdb,1,basis)
107
 
c      call schwarz_init(geom,basis)
 
121
      call int_init(rtdb,1,basis)
 
122
      call schwarz_init(geom,basis)
108
123
      call hnd_giao_init(basis,1)
109
124
      call scf_get_fock_param(rtdb,tol2e)
110
 
c
111
125
      status = rtdb_parallel(.true.)
112
 
 
113
126
c     Get Unperturbed MO vectors and eigenvalues
114
127
c     First allocate some memory for occupation numbers and eigenvalues
115
 
c
116
128
      if (.not. bas_numbf(basis,nbf)) call
117
129
     &    errquit('giao_b1: could not get nbf',0, BASIS_ERR)
118
130
      if (.not. ma_push_get(mt_dbl,2*nbf,'occ num',l_occ,k_occ)) call
119
131
     &    errquit('giao_b1: ma_push_get failed k_occ',0,MA_ERR)
120
132
      if (.not. ma_push_get(mt_dbl,2*nbf,'eigenval',l_eval,k_eval)) call
121
133
     &    errquit('giao_b1: ma_push_get failed k_eval',0,MA_ERR)
122
 
      call hnd_prp_vec_read(rtdb,geom,basis,nbf,nclosed,nopen,
123
 
     &                      nvirt,scftyp,vectors,dbl_mb(k_occ),
124
 
     &                      dbl_mb(k_eval),nmo)
 
134
 
 
135
      call hnd_prp_vec_read(rtdb,geom,basis,     ! in : handles
 
136
     &                      nbf,                 ! out: nr basis functions
 
137
     &                      nclosed,nopen,nvirt, ! out: occupation numbers
 
138
     &                      scftyp,              ! out: type calc
 
139
     &                      vectors,             ! out: MO vectors
 
140
     &                      dbl_mb(k_occ),       ! out: occupations
 
141
     &                      dbl_mb(k_eval),      ! out: DFT energies
 
142
     &                      nmo)                 ! out: nr MOs
 
143
 
 
144
      call get_nocc(rtdb,   ! in : rtdb handle
 
145
     &              nocc,   ! out: nr occupations
 
146
     &              npol,   ! out: nr of polarization
 
147
     &              nclosed,! in : nr closed shells
 
148
     &              nopen,  ! in : nr open shells
 
149
     &              nvirt,  ! in : nr virtual MOs
 
150
     &              scftyp, ! in : string = UHF or RHF
 
151
     &              ntot)   ! out: sum_{i,npol} nocc(i)*nvirt(i)
 
152
 
 
153
      if (ga_nodeid().eq.0) then
 
154
        write(*,10) nocc(1)   ,nocc(2),
 
155
     &              nopen(1)  ,nopen(2),
 
156
     &              nclosed(1),nclosed(2),
 
157
     &              nvirt(1)  ,nvirt(2),scftyp,ntot
 
158
 10    format('giao_b1_mov: nocc =(',i3,',',i3,') ',
 
159
     &        'nopen=(',i3,',',i3,') ',
 
160
     &        'nclos=(',i3,',',i3,') ',
 
161
     &        'nvirt=(',i3,',',i3,') ',
 
162
     &        'scftyp=',a,' ntot=',i3)
 
163
      endif
125
164
c
126
165
c     Get Unperturbed Density Matrix
127
 
c
128
 
      call hnd_prp_get_dens(rtdb,geom,basis,g_dens,ndens,scftyp,
129
 
     &                      nclosed,nopen,nvirt)
 
166
 
 
167
      call hnd_prp_get_dens(rtdb,geom,basis,     ! in : handles
 
168
     &                      g_dens,ndens,        ! out: electron density
 
169
     &                      scftyp,              ! in : type calc
 
170
     &                      nclosed,nopen,nvirt) ! in : occupations numbers
 
171
 
 
172
c      if (ga_nodeid().eq.0) then
 
173
c        write(*,21) npol,nocc(1)   ,nocc(2),
 
174
c     &              nopen(1)  ,nopen(2),
 
175
c     &              nclosed(1),nclosed(2),
 
176
c     &              nvirt(1)  ,nvirt(2),scftyp,ntot
 
177
c 21    format('npol=',i3,' nocc =(',i3,',',i3,') ',
 
178
c     &        'nopen=(',i3,',',i3,') ',
 
179
c     &        'nclos=(',i3,',',i3,') ',
 
180
c     &        'nvirt=(',i3,',',i3,') ',
 
181
c     &        'scftyp=',a,' ntot=',i3)
 
182
c      endif
130
183
 
131
184
      if (debug) write (luout,*) 'unpertubed MOs and Pmat assembled'
132
185
c
133
 
c     Error exit if scftyp equals UHF (= ROHF)
134
 
c
135
 
      if (scftyp.eq.'UHF') then
136
 
          if (ga_nodeid().eq.0) write(luout,7000)
137
 
          call errquit('giao_b1: incompatible SCF type for Response',
138
 
     &       0,INPUT_ERR)
139
 
      endif
140
 
c
141
186
c     Create U matrix of dimension (nbf,nmo,3) and zero
142
187
c     Use ahi for dimension and ahi array for chunking/blocking
143
 
c
144
 
      alo(1) = nbf
145
 
      alo(2) = -1
146
 
      alo(3) = -1
147
 
      ahi(1) = nbf
148
 
      ahi(2) = nclosed(1)
149
 
      ahi(3) = 3
150
 
      if (.not.nga_create(MT_DBL,3,ahi,'U matrix',alo,g_u)) call 
 
188
      do ispin=1,npol
 
189
       alo(1) = nbf
 
190
       alo(2) = -1
 
191
       alo(3) = -1
 
192
       ahi(1) = nbf
 
193
       ahi(2) = nocc(ispin)
 
194
       ahi(3) = ndir
 
195
       if (.not.nga_create(MT_DBL,3,ahi,'U matrix',alo,g_u(ispin))) call 
151
196
     &    errquit('giao_b1: nga_create failed g_u',0,GA_ERR)
152
 
      call ga_zero(g_u)
153
 
c
 
197
       call ga_zero(g_u(ispin))
 
198
      enddo ! end-loop-ispin
 
199
c 
154
200
c     Construction of right-hand side CPHF
155
201
c     Create CPHF array of proper dimension : (nocc*nvirt,3)
156
 
c
157
 
      if(.not.ga_create(MT_DBL,nclosed(1)*nvirt(1),3,'RHS',-1,-1,g_rhs))
 
202
      ndata=2 ! 1st subspace corresponds to g_b, 
 
203
c             ! 2nd subspace corresponds to sol (if exists)
 
204
        if (.not. rtdb_put(rtdb,'cphf2-aores:ndata', 
 
205
     &          mt_int, 1,ndata)) call
 
206
     $     errquit('fiao_b1: failed to write skew ', 0, RTDB_ERR)
 
207
      if(.not.ga_create(MT_DBL,ntot,ndata*ndir,
 
208
     &   'RHS',-1,-1,g_rhs))
158
209
     &   call errquit('giao_b1: ga_create failed g_rhs',0,GA_ERR)
159
210
      call ga_zero(g_rhs)
160
 
c
161
 
c     NGA dimension arrays for copying will be the same every time
162
 
c     Also third NGA dimension for any of the three dimensional
163
 
c     arrays will be the same everytime (running from 1 to 3)
164
 
c     So, lets define them once and for all in blo and bhi
165
 
c
166
 
      blo(1) = 1
167
 
      bhi(1) = nclosed(1)*nvirt(1)
168
 
      blo(2) = 1
169
 
      bhi(2) = 3
170
 
c    
171
 
c     Get S10 in GA and transform to MO set (virt,occ)
172
 
c
173
 
      alo(1) = nbf
174
 
      alo(2) = -1
175
 
      alo(3) = -1
176
 
      ahi(1) = nbf
177
 
      ahi(2) = nbf
178
 
      ahi(3) = 3
179
 
      if (.not.nga_create(MT_DBL,3,ahi,'s10 matrix',alo,g_s10)) call 
180
 
     &    errquit('giao_b1: nga_create failed g_s01',0,GA_ERR)
181
 
      call ga_zero(g_s10)
182
 
      call int_giao_1ega(basis,basis,g_s10,'s10',origin,
183
 
     &                   nat,oskel)
184
 
      call giao_aotomo(g_s10,vectors(1),nclosed,nvirt,1,3,nbf)
185
 
 
186
 
      if (debug) write (luout,*) 'S10 done'
187
 
 
188
 
c     while we are calculating integrals, let's also determine
189
 
c     the 'half' overlap derivative, used later in the calling
190
 
c     routine
191
 
 
192
 
      call ga_zero(g_sket1)
193
 
      call int_giao_1ega(basis,basis,g_sket1,'srxRb',origin,
194
 
     &   nat,oskel)
195
 
 
196
 
      if (debug) write (luout,*) 'S1-ket done'
197
 
 
198
 
c     g_sket1 will not be used further here. It is one of the 
199
 
c     output results of this routine.
200
 
c     Broceed with the computation of the B-field perturbed
201
 
c     MO coefficients.
202
 
 
203
 
c
204
 
c     ga_rhs(a,i) = ga_rhs(a,i) - e(i) * S10(a,i)
205
 
c     Scale (occ,virt) block g_s10 with - (minus) eigenvalues 
206
 
c
207
 
      alo(1) = nclosed(1)+1
208
 
      ahi(1) = nmo
209
 
      alo(3) = 1
210
 
      ahi(3) = 3
211
 
      do iocc = 1, nclosed(1)
212
 
         alo(2) = iocc
213
 
         ahi(2) = iocc
214
 
         call nga_scale_patch(g_s10,alo,ahi,-dbl_mb(k_eval+iocc-1)) 
215
 
      enddo
216
 
      if (.not.ma_pop_stack(l_eval)) call
217
 
     &    errquit('giao_b1: ma_pop_stack failed k_eval',0,MA_ERR)
218
 
      if (.not.ma_pop_stack(l_occ)) call
219
 
     &    errquit('giao_b1: ma_pop_stack failed k_occ',0,MA_ERR)
220
 
c
221
 
c     Copy to ga_rhs 
222
 
c     alo(1) and ahi(1) the same as before
223
 
c
224
 
      alo(2) = 1
225
 
      ahi(2) = nclosed(1)
226
 
      call nga_copy_patch('n',g_s10,alo,ahi,g_rhs,blo,bhi)
227
 
 
228
 
      if (debug) write (luout,*) 'S10 to Umat done'
229
 
c
230
 
c     Construct occ-occ part of the three U matrices
231
 
c     Occ-occ blocks for each field direction are defined as -1/2 S10
232
 
c     Scale (occ,occ) block g_s10 with -1/2 and add to g_u
233
 
c
234
 
c     alo(2) and ahi(2) will stay as 1 and nclosed(1) for a while
235
 
c
236
 
      alo(1) = 1
237
 
      ahi(1) = nclosed(1)
238
 
      call nga_scale_patch(g_s10,alo,ahi,-0.5d0)
239
 
      call nga_copy_patch('n',g_s10,alo,ahi,g_u,alo,ahi)
240
 
 
241
 
      if (debug) write (luout,*) 'S10 in occ-occ done'
 
211
 
 
212
c      if (ga_nodeid().eq.0)
 
213
c     &  write(*,*) 'FA-BEF update_rhs_eS10'
 
214
 
 
215
c      if (ga_nodeid().eq.0) then
 
216
c        write(*,70) nocc(1)   ,nocc(2),
 
217
c     &              nvirt(1)  ,nvirt(2),npol,nbf,nmo
 
218
c 70    format('BEF-update_rhs_eS10 nocc =(',i5,',',i5,') ',
 
219
c     &        'nvirt=(',i5,',',i5,') ',
 
220
c     &        '(npol,nbf,nmo)=(',i3,',',i3,',',i3,')')
 
221
c      endif
 
222
 
 
223
      call update_rhs_eS10(
 
224
     &            g_rhs,         !in/out:
 
225
     &            g_u,           !out:
 
226
     &            g_sket1,       !out:
 
227
     &            dbl_mb(k_eval),!in : energy values
 
228
     &            vectors,       !in : MO vectors
 
229
     &            nocc,          !in : nr.   occupied MOs
 
230
     &            nvirt,         !in : nr. unoccupied MOs
 
231
     &            npol,          !in : nr. polarizations
 
232
     &            nbf,           !in : nr. basis functions
 
233
     &            nmo,           !in : nr. MOs
 
234
     &            basis,         !in : basis handle
 
235
     &            debug)         !in : logical var for debugging
 
236
 
 
237
      if (debug) then
 
238
       if (ga_nodeid().eq.0)
 
239
     &  write(*,*) 'FA-AFT update_rhs_eS10'
 
240
       if (ga_nodeid().eq.0)
 
241
     &  write(*,*) '---- g_rhs-AFT-eS10-------- START'
 
242
        call ga_print(g_rhs)
 
243
       if (ga_nodeid().eq.0)
 
244
     &  write(*,*) '---- g_rhs-AFT-eS10--------  END'
 
245
       if (ga_nodeid().eq.0)
 
246
     &  write(*,*) '---- g_u-AFT-eS10-------- START'
 
247
       do ispin=1,npol
 
248
         call ga_print(g_u(ispin))
 
249
       enddo ! end-loop-ispin
 
250
       if (ga_nodeid().eq.0)
 
251
     &  write(*,*) '---- g_u-AFT-eS10--------  END'
 
252
       if (ga_nodeid().eq.0)
 
253
     &  write(*,*) '---- g_sket1-AFT-eS10-------- START'
 
254
        call ga_print(g_sket1)
 
255
       if (ga_nodeid().eq.0)
 
256
     &  write(*,*) '---- g_sket1-AFT-eS10--------  END'
 
257
      endif ! end-if-debug
 
258
 
242
259
c
243
260
c     We also need the occupied-occupied contribution of g_u
244
261
c     contributing to the first order density matrix. As this block does
245
262
c     not change during the CPHF we can calculate it once and subtract
246
263
c     it from the RHS. We will reuse g_s10 as scratch space.
247
 
c
248
 
      call ga_zero(g_s10)
249
 
      clo(1) = 3
250
 
      clo(2) = nbf
251
 
      clo(3) = nbf
252
 
      chi(1) = 1  
253
 
      chi(2) = -1 
254
 
      chi(3) = -1
255
 
      if (.not.nga_create(MT_DBL,3,clo,'Fock matrix',chi,g_fock)) call 
256
 
     &    errquit('giao_b1: nga_create failed g_fock',0,GA_ERR)
257
 
      if (.not.nga_create(MT_DBL,3,clo,'D10 matrix',chi,g_d1)) call 
258
 
     &    errquit('giao_b1: nga_create failed g_d1',0,GA_ERR)
259
 
      call ga_zero(g_fock)
260
 
      call ga_zero(g_d1)
261
 
      alo(1) = 1
262
 
      alo(2) = 1
263
 
      blo(1) = 1
264
 
      blo(2) = 1
265
 
      bhi(1) = nbf
266
 
      clo(2) = 1
267
 
      clo(3) = 1
268
 
      chi(2) = nbf
269
 
      chi(3) = nbf
270
 
      dlo(1) = 1
271
 
      dlo(2) = 1
272
 
      dhi(1) = nbf
273
 
      dhi(2) = nclosed(1)
274
 
c
275
 
c     Create "perturbed density matrix" for closed-closed g_u block
276
 
c
277
 
      do ifld = 1, 3
278
 
         jfac(ifld) = 0.0d0
279
 
         kfac(ifld) = -1.0d0*xfac
280
 
         alo(3) = ifld
281
 
         ahi(3) = ifld
282
 
         clo(1) = ifld
283
 
         chi(1) = ifld
284
 
         dlo(3) = ifld
285
 
         dhi(3) = ifld
286
 
         ahi(1) = nmo
287
 
         ahi(2) = nclosed(1)
288
 
         bhi(2) = nmo 
289
 
         call nga_matmul_patch('n','n',1.0d0,0.0d0,vectors(1),blo,bhi,  
290
 
     &                         g_u,alo,ahi,g_s10,dlo,dhi)  
291
 
         ahi(2) = nbf
292
 
         ahi(1) = nclosed(1)
293
 
         bhi(2) = nclosed(1)
294
 
c
295
 
c     Minus sign as we subtract it from the RHS as we do not include 
296
 
c     it in the LHS
297
 
c
298
 
         call nga_matmul_patch('n','t',-1.0d0,0.0d0,vectors(1),blo,bhi,
299
 
     &                         g_s10,alo,ahi,g_d1,clo,chi)  
300
 
      enddo
301
 
c
302
 
c     Build "first order fock matrix"
303
 
c
304
 
      if (use_theory.eq.'dft') then
305
 
         if(.not. rtdb_put(rtdb,'bgj:xc_active', MT_LOG, 1, .true.))
306
 
     $     call errquit('hess_cphf: rtdb_put of xc_active failed',0,
307
 
     &       RTDB_ERR)
308
 
         if(.not. rtdb_put(rtdb,'fock_xc:calc_type', MT_INT, 1, 2))
309
 
     $     call errquit('hess_cphf: rtdb_put of calc_type failed',0,
310
 
     &       RTDB_ERR)
311
 
         if(.not. rtdb_put(rtdb,'fock_j:derfit', MT_LOG, 1, .false.))
312
 
     $     call errquit('hess_cphf: rtdb_put of j_derfit failed',0,
313
 
     &       RTDB_ERR)
314
 
      endif
315
 
c      call shell_fock_build(geom, basis, 0, 3,
316
 
c     $     jfac, kfac,tol2e, g_d1, g_fock,.false.)
317
 
c      if(use_theory.eq.'dft') then
318
 
c         if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, 0))
319
 
c     $      call errquit('giaox: rtdb_put failed',0,RTDB_ERR)
 
264
 
 
265
      call get_d1_giao_b1(g_d1,   ! out:
 
266
     &                    g_u,    ! in :
 
267
     &                    vectors,! in : MO vectors
 
268
     &                    nocc,   ! in : nr. occ shells
 
269
     &                    npol,   ! in : nr. polarizations
 
270
     &                    nbf,    ! in : nr. basis functions
 
271
     &                    nmo,    ! in : nr. MOs
 
272
     &                    debug)  ! in : =.true. for debugging
 
273
 
 
274
      if (debug) then
 
275
       if (ga_nodeid().eq.0)
 
276
     &  write(*,*) '------- g_d1-aft-get_d1 ---- START'
 
277
        call ga_print(g_d1)
 
278
       if (ga_nodeid().eq.0)
 
279
     &  write(*,*) '------- g_d1-aft-get_d1 ---- END'
 
280
      endif ! end-if-debug
 
281
c      if (ga_nodeid().eq.0)
 
282
c     & write(*,*) 'FA-BEF update_rhs_shfock'
 
283
c      if (ga_nodeid().eq.0) then
 
284
c        write(*,71) nocc(1)   ,nocc(2),
 
285
c     &              nvirt(1)  ,nvirt(2),npol,nbf,nmo
 
286
c 71    format('BEF-update_rhs_shfock =(',i5,',',i5,') ',
 
287
c     &        'nvirt=(',i5,',',i5,') ',
 
288
c     &        '(npol,nbf,nmo)=(',i3,',',i3,',',i3,')')
320
289
c      endif
321
290
 
322
 
c from hnd_giaox.F]
323
 
c
324
 
c     Note: Just the exchange: jfac = 0.d0 (see above)
325
 
c
326
 
      if (.not.cam_exch) then
327
 
         call shell_fock_build(geom, basis, 0, 3,
328
 
     $     jfac, kfac, tol2e, g_d1, g_fock, .false.)
329
 
      else
330
 
         call shell_fock_build_cam(geom, basis, 0, 3,
331
 
     $     jfac, kfac, tol2e, g_d1, g_fock, .false.)
332
 
      end if
333
 
c
334
 
      if(use_theory.eq.'dft') then
335
 
         if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, 0))
336
 
     $      call errquit('giaox: rtdb_put failed',0,RTDB_ERR)
337
 
      endif
338
 
 
339
 
 
340
 
c
341
 
c     Transform to the occ-virt MO basis and add to RHS
342
 
c
343
 
      call ga_zero(g_s10)
344
 
      alo(1) = 1
345
 
      ahi(1) = nbf
346
 
      alo(2) = 1
347
 
      ahi(2) = nclosed(1)
348
 
      do ifld = 1, 3
349
 
         alo(3) = ifld
350
 
         ahi(3) = ifld
351
 
         clo(1) = ifld
352
 
         chi(1) = ifld
353
 
         call nga_matmul_patch('n','n',2.0d0,0.0d0,
354
 
     $        g_fock,clo,chi,
355
 
     $        vectors(1),alo,ahi,
356
 
     $        g_s10,alo,ahi)
357
 
      enddo
358
 
      call ga_zero(g_fock)
359
 
      clo(2) = nclosed(1)+1
360
 
      clo(3) = 1
361
 
      chi(2) = nmo
362
 
      chi(3) = nclosed(1)
363
 
      do ifld = 1, 3
364
 
         blo(1) = nclosed(1)+1
365
 
         blo(2) = 1
366
 
         bhi(1) = nmo
367
 
         bhi(2) = nbf
368
 
         alo(3) = ifld
369
 
         ahi(3) = ifld
370
 
         clo(1) = ifld
371
 
         chi(1) = ifld
372
 
         call nga_matmul_patch('t','n',1.0d0,0.0d0,
373
 
     $        vectors(1), blo,bhi,
374
 
     $        g_s10, alo,ahi,
375
 
     $        g_fock, clo,chi )
376
 
         blo(1) = 1
377
 
         bhi(1) = nclosed(1)*nvirt(1)
378
 
         blo(2) = ifld
379
 
         bhi(2) = ifld
380
 
         call nga_add_patch(1.0d0,g_rhs,blo,bhi,1.0d0,g_fock,clo,chi,
381
 
     &                      g_rhs,blo,bhi)
382
 
 
383
 
!      call nga_print_patch(g_fock,clo,chi,1)
384
 
 
385
 
      enddo
386
 
c
387
 
c     Cleanup of g_d1 and g_fock, not needed for now
388
 
c
 
291
      call update_rhs_shfock(
 
292
     &                    g_rhs,  ! in/out: RHS used for cphf2/3
 
293
     &                    g_d1,   ! in    :
 
294
     &                    vectors,! in    : MO vectors
 
295
     &                    rtdb,   ! in    : rtdb  handle
 
296
     &                    geom,   ! in    : geom  handle
 
297
     &                    basis,  ! in    : basis handle 
 
298
     &                    jfac,   ! in    : exch factors
 
299
     &                    kfac,   ! in    : exch factors
 
300
     &                    tol2e,  ! in    : tolerance coeff
 
301
     &                    nocc,   ! in    : nr. occ  shells
 
302
     &                    nvirt,  ! in    : nr. virt shells
 
303
     &                    npol,   ! in    : nr. polarizations
 
304
     &                    nbf,    ! in    : nr. basis functions
 
305
     &                    nmo,    ! in    : nr. MOs
 
306
     &                    debug)  ! in    : =.true. for debugging
 
307
 
 
308
      if (debug) then
 
309
       if (ga_nodeid().eq.0)
 
310
     &  write(*,*) 'FA-AFT update_rhs_shfock'
 
311
       if (ga_nodeid().eq.0)
 
312
     &  write(*,*) '---- g_rhs-AFT-shfock-------- START'
 
313
        call ga_print(g_rhs)
 
314
       if (ga_nodeid().eq.0)
 
315
     &  write(*,*) '---- g_rhs-AFT-shfock--------  END'
 
316
      endif ! end-if-debug
 
317
 
389
318
      if (.not.ga_destroy(g_d1)) call 
390
319
     &    errquit('giao_b1: ga_destroy failed g_d1',0,GA_ERR)
391
 
      if (.not.ga_destroy(g_fock)) call 
392
 
     &    errquit('giao_b1: ga_destroy failed g_fock',0,GA_ERR)
393
320
 
394
321
      if (debug) write (luout,*) 'S10 Fockop done'
395
 
c
396
 
c     Get H10 in GA, reusing g_s10 array
397
 
c
398
 
      call ga_zero(g_s10)
399
 
      call int_giao_1ega(basis,basis,g_s10,'l10',origin,
400
 
     &                   nat,oskel)
401
 
      call int_giao_1ega(basis,basis,g_s10,'tv10',origin,
402
 
     &                   nat,oskel)
403
 
c
404
 
c     Get external and cosmo bq contribution
405
 
c
406
 
      nbq = 0
407
 
      nextbq = 0
408
 
      ncosbq = 0
409
 
      if(geom_extbq_on()) nextbq = geom_extbq_ncenter()
410
 
      nbq = nextbq ! external bq's
411
 
c
412
 
      if (rtdb_get(rtdb,'cosmo:nefc',mt_int,1,ncosbq))
413
 
     &    nbq = ncosbq ! cosmo bq's
414
 
c
415
 
      if (nextbq.gt.0.and.ncosbq.gt.0)
416
 
     &    nbq = nextbq + ncosbq  ! tally up cosmo and external bqs
417
 
c
418
 
c     if (ga_nodeid().eq.0) write(6,*) "nbq: ", nbq
419
 
      if (nbq.gt.0) then
420
 
        call int_giao_1ega(basis,basis,g_s10,'bq10',origin,
421
 
     &                   nat,oskel)
422
 
      end if
423
 
c
424
 
c     ga_rhs(a,i) = ga_rhs(a,i) + H10(a,i)
425
 
c     Transform H10 to MO and add to g_rhs
426
 
c
427
 
      call giao_aotomo(g_s10,vectors(1),nclosed,nvirt,1,3,nbf)
428
 
      alo(1) = nclosed(1)+1
429
 
      ahi(1) = nmo
430
 
      alo(2) = 1
431
 
      ahi(2) = nclosed(1)
432
 
      alo(3) = 1
433
 
      ahi(3) = 3
 
322
c      if (ga_nodeid().eq.0)
 
323
c     & write(*,*) 'FA-BEF update_rhs_threeAOints'
 
324
 
 
325
      call update_rhs_threeAOints(
 
326
     &                    g_rhs,  ! in/out: RHS used for cphf2/3
 
327
     &                    vectors,! in    : MO vectors
 
328
     &                    rtdb,   ! in    : rtdb  handle 
 
329
     &                    basis,  ! in    : basis handle 
 
330
     &                    nocc,   ! in    : nr occ  shells
 
331
     &                    nvirt,  ! in    : nr virt shells
 
332
     &                    npol,   ! in    : nr. polarizations
 
333
     &                    nbf,    ! in    : nr. basis functions
 
334
     &                    nmo,    ! in    : nr. MOs
 
335
     &                    debug)  ! in    : logical for debugging
 
336
 
 
337
      if (debug) then
 
338
       if (ga_nodeid().eq.0)
 
339
     &  write(*,*) 'FA-AFT update_rhs_threeAOints'
 
340
       if (ga_nodeid().eq.0)
 
341
     &  write(*,*) '---- g_rhs-AFT-threeAOints-------- START'
 
342
        call ga_print(g_rhs)
 
343
       if (ga_nodeid().eq.0)
 
344
     &  write(*,*) '---- g_rhs-AFT-threeAOints--------  END'
 
345
      endif ! end-if-debug
 
346
 
 
347
      call update_rhs_fock2e(
 
348
     &                    g_rhs,  ! in/out: 
 
349
     &                    vectors,! in : MO vectors
 
350
     &                    rtdb,   ! in : rtdb  handle
 
351
     &                    basis,  ! in : basis handle
 
352
     &                    geom,   ! in : geom  handle
 
353
     &                    g_dens, ! in : e-density
 
354
     &                    nocc,   ! in : nr. occ  shells
 
355
     &                    nvirt,  ! in : nr. virt shells
 
356
     &                    npol,   ! in : nr. polarizations
 
357
     &                    nbf,    ! in : nr. basis functions
 
358
     &                    nmo,    ! in : nr. MOs   
 
359
     &                    xfac,   ! in : exchange factor
 
360
     &                    tol2e,  ! in : tolerance coeff.
 
361
     &                    debug)  ! in : logical for debugging
 
362
 
 
363
      call schwarz_tidy()
 
364
      call int_terminate()
 
365
 
434
366
      blo(1) = 1
435
 
      bhi(1) = nclosed(1)*nvirt(1)
 
367
      bhi(1) = ntot
436
368
      blo(2) = 1
437
 
      bhi(2) = 3
438
 
      call nga_add_patch(1.0d0,g_rhs,blo,bhi,1.0d0,g_s10,alo,ahi,
439
 
     &                   g_rhs,blo,bhi)
440
 
c
441
 
c     Cleanup g_s10 as we do not need it right now
442
 
c
443
 
      if (.not.ga_destroy(g_s10)) call 
444
 
     &    errquit('giao_b1: ga_destroy failed g_s10',0,GA_ERR)
445
 
 
446
 
      if (debug) write (luout,*) 'H10 to Umat done'
447
 
c
448
 
c     Remaining term is Perturbed (GIAO) two-electron term times
449
 
c     Unperturbed density Calculate Sum(r,s) D0(r,s) * G10(m,n,r,s) in
450
 
c     AO basis
451
 
c
452
 
      alo(1) = -1 
453
 
      alo(2) = -1
454
 
      alo(3) = 1
455
 
      ahi(1) = nbf
456
 
      ahi(2) = nbf
457
 
      ahi(3) = 3
458
 
      if (.not.nga_create(MT_DBL,3,ahi,'Fock matrix',alo,g_fock)) call 
459
 
     &    errquit('giao_b1: nga_create failed g_fock',0,GA_ERR)
460
 
      call ga_zero(g_fock)
461
 
      if(use_theory.eq.'dft') then
462
 
         ifld = 4
463
 
         if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, ifld))
464
 
     $      call errquit('giao_b1: rtdb_put failed',0,RTDB_ERR)
465
 
      endif
466
 
      if (debug) write (luout,*) 'calling new_giao_2e'
467
 
c     call new_giao_2e(geom, basis, nbf, tol2e, g_dens, g_fock, xfac)
468
 
      call new_giao_2e(geom, basis, nbf, tol2e, g_dens, g_fock, xfac,1) ! FA-restricted calc
469
 
      if (debug) write (luout,*) 'done new_giao_2e'
470
 
      if(use_theory.eq.'dft') then
471
 
         ifld = 0
472
 
         if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, ifld))
473
 
     $      call errquit('giao_b1: rtdb_put failed',0,RTDB_ERR)
474
 
         if(.not. rtdb_put(rtdb,'bgj:xc_active', MT_LOG, 1, .false.))
475
 
     $       call errquit('giao_b1: rtdb_put of xc_active failed',0,
476
 
     &       RTDB_ERR)
477
 
      endif
478
 
c
479
 
c     Transform to MO basis and add to right-hand-side
480
 
c
481
 
      call giao_aotomo(g_fock,vectors(1),nclosed,nvirt,1,3,nbf)
482
 
      alo(1) = nclosed(1)+1
483
 
      ahi(1) = nmo
484
 
      alo(2) = 1
485
 
      ahi(2) = nclosed(1)
486
 
      alo(3) = 1
487
 
      ahi(3) = 3
488
 
      call nga_add_patch(1.0d0,g_rhs,blo,bhi,1.0d0,g_fock,alo,ahi,
489
 
     &                   g_rhs,blo,bhi)
490
 
      if (.not.ga_destroy(g_fock)) call 
491
 
     &    errquit('giao_b1: ga_destroy failed g_fock',0,GA_ERR)
492
 
      call nga_scale_patch(g_rhs,blo,bhi,-4.0d0)
493
 
 
494
 
      if (debug) write (luout,*) 'D00 GIAO Fock terms done'
 
369
      bhi(2) = ndir  
 
370
      if      (npol.eq.1) then
 
371
        call nga_scale_patch(g_rhs,blo,bhi,-4.0d0)
 
372
      else if (npol.eq.2) then
 
373
        call nga_scale_patch(g_rhs,blo,bhi,-2.0d0)
 
374
      endif  
 
375
 
 
376
      if (debug) then
 
377
       if (ga_nodeid().eq.0)
 
378
     &  write(*,*) 'FA-AFT update_rhs_fock2e'
 
379
       if (ga_nodeid().eq.0)
 
380
     &  write(*,*) '---- g_rhs-AFT-fock2e-------- START'
 
381
        call ga_print(g_rhs)
 
382
       if (ga_nodeid().eq.0)
 
383
     &  write(*,*) '---- g_rhs-AFT-fock2e--------  END'
 
384
      endif ! end-if-debug
 
385
 
 
386
       call util_file_name(lbl_cphfaoresp,
 
387
     &                     .false.,.false.,aorespfilename)
 
388
      read_grhs_giaob1=0 
 
389
      if (.not. dft_CPHF1_read( ! file exists and read g_rhs guess
 
390
     &           aorespfilename,! in: filename
 
391
     &           npol,          ! in: nr polarization
 
392
     &           nocc,          ! in: nr occupied MOs
 
393
     &           nvirt,         ! in: nr virtual  MOs
 
394
     &           1,             ! in: nr. components
 
395
     &           g_rhs,         ! in: (ntot,3)       GA matrix
 
396
     &           g_rhs_im,      ! in: dummy
 
397
     &           .false.))      ! in: =T if (RE,IM) =F if RE
 
398
     & then
 
399
         read_grhs_giaob1=1
 
400
       else
 
401
        if (.not. rtdb_put(rtdb,'cphf2-aores:guess', 
 
402
     &          mt_log, 1,.true.)) call
 
403
     $     errquit('giao_b1: failed to write skew ', 0, RTDB_ERR)
 
404
      endif
 
405
 
 
406
      if (debug) then
 
407
       if (ga_nodeid().eq.0)
 
408
     &  write(*,*) '---- g_rhs-AFT-readfile-------- START'
 
409
        call ga_print(g_rhs)
 
410
       if (ga_nodeid().eq.0)
 
411
     &  write(*,*) '---- g_rhs-AFT-readfile--------  END'
 
412
      endif ! end-if-debug
 
413
  
 
414
c       if (ga_nodeid().eq.0)
 
415
c     &  write(*,*) 'COMPUTE cphf giao_b1 data ...'   
495
416
c
496
417
c     Write ga_rhs to disk 
497
 
c
498
418
      call util_file_name('cphf_rhs',.true.,.true.,cphf_rhs)
499
419
      call util_file_name('cphf_sol',.true.,.true.,cphf_sol)
500
420
      if(.not.file_write_ga(cphf_rhs,g_rhs)) call errquit
501
421
     $  ('giao_b1: could not write cphf_rhs',0, DISK_ERR)
502
422
c
503
 
      call schwarz_tidy()
504
 
      call int_terminate()
505
 
c
506
423
c     Call the CPHF routine
507
424
c     
508
425
c     We do need to tell the CPHF that the density is skew symmetric.
509
426
c     Done via rtdb, put cphf:skew .false. on rtdb and later remove it.
510
 
c
511
427
      if (.not. rtdb_put(rtdb, 'cphf:skew', mt_log, 1,.false.)) call
512
428
     $   errquit('giao_b1: failed to write skew ', 0, RTDB_ERR)
513
 
 
514
429
      if (debug) write (luout,*) 'calling cphf'
515
 
c
516
430
      if (.not.cphf2(rtdb)) call errquit
517
431
     $  ('giao_b1: failure in cphf ',0, RTDB_ERR)
518
 
c
519
432
      if (.not. rtdb_delete(rtdb, 'cphf:skew')) call
520
433
     $   errquit('giao_b1: rtdb_delete failed ', 0, RTDB_ERR)
521
 
 
522
434
      if (debug) write (luout,*) 'cphf done'
523
435
c
524
436
c     Occ-virt blocks are the solution pieces of the CPHF
525
437
c     Read solution vector from disk and put solutions in U matrices
526
 
c
527
438
      call ga_zero(g_rhs)
528
439
      if(.not.file_read_ga(cphf_sol,g_rhs)) call errquit
529
 
     $  ('giao_b1: could not read cphf_rhs',0, DISK_ERR)      
530
 
      call nga_copy_patch('n',g_rhs,blo,bhi,g_u,alo,ahi)
531
 
c
532
 
      if (.not.ga_destroy(g_rhs)) call 
533
 
     &    errquit('giao_b1: ga_destroy failed g_rhs',0,GA_ERR)
534
 
c
535
 
c     From U matrices, generate the perturbed density matrices D1x,y,z
536
 
c     C1 = C0 * U10
537
 
c     D1 = 2[(C1*C0+) - (C0*C1+)]
538
 
c
539
 
      alo(1) = nbf
540
 
      alo(2) = -1
541
 
      alo(3) = -1
542
 
      ahi(1) = nbf
543
 
      ahi(2) = nbf
544
 
      ahi(3) = 3
545
 
      if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1)) call 
546
 
     &    errquit('giao_b1: nga_create failed g_d1',0,GA_ERR)
547
 
c
548
 
      alo(1) = 1
549
 
      alo(2) = 1
550
 
      blo(1) = 1
551
 
      blo(2) = 1
552
 
      clo(1) = 1
553
 
      chi(1) = nbf
554
 
      clo(2) = 1
555
 
      chi(2) = nbf
556
 
      dlo(1) = 1
557
 
      dlo(2) = 1
558
 
      dhi(1) = nbf
559
 
      dhi(2) = nclosed(1)
560
 
      do ifld = 1, 3
561
 
         alo(3) = ifld
562
 
         ahi(3) = ifld
563
 
         blo(3) = ifld
564
 
         bhi(3) = ifld
565
 
         clo(3) = ifld
566
 
         chi(3) = ifld
567
 
         dlo(3) = ifld
568
 
         dhi(3) = ifld
569
 
         bhi(1) = nbf
570
 
         bhi(2) = nmo 
571
 
         ahi(1) = nmo
572
 
         ahi(2) = nclosed(1)
573
 
c
574
 
c     Make C1
575
 
c
576
 
         call nga_matmul_patch('n','n',1.0d0,0.0d0,vectors(1),blo,bhi,  
577
 
     &                         g_u,alo,ahi,g_d1,dlo,dhi)  
578
 
         call nga_copy_patch('n',g_d1,dlo,dhi,g_vecB1,dlo,dhi)
579
 
 
580
 
c        This patch of g_vecB1 now has the perturbed MO
581
 
c        coefficients. let's print them for debug purposes:
582
 
 
583
 
         if (.false.) then ! change .false. to debug to print them
584
 
           write (luout,*) 'giao_b1: perturbed C, direction ',ifld
585
 
           call nga_print_patch(g_vecB1,dlo,dhi,1)
586
 
         end if
587
 
 
588
 
       enddo ! ifld
589
 
 
590
 
      if (.not.ga_destroy(g_u)) call 
591
 
     &    errquit('giao_b1: ga_destroy failed g_u',0,GA_ERR)
592
 
      if (.not.ga_destroy(g_d1)) call 
593
 
     &    errquit('giao_b1: ga_destroy failed g_d1',0,GA_ERR)
594
 
      if (.not.ga_destroy(vectors(1))) call 
595
 
     &   errquit('giao_b1: ga_destroy failed vectors',0,GA_ERR)
596
 
      if (.not.ga_destroy(g_dens(1))) call 
 
440
     $  ('giao_b1: could not read cphf_rhs',0, DISK_ERR)  
 
441
 
 
442
      if (debug) then
 
443
       if (ga_nodeid().eq.0)
 
444
     &  write(*,*) '---- g_rhs-BEF-write2file------- START'
 
445
        call ga_print(g_rhs)
 
446
       if (ga_nodeid().eq.0)
 
447
     &  write(*,*) '---- g_rhs-BEF-write2file--------  END'
 
448
      endif ! end-if-debug
 
449
 
 
450
       call util_file_name(lbl_cphfaoresp,
 
451
     &                     .false.,.false.,aorespfilename)
 
452
 
 
453
       status=dft_CPHF1_write(
 
454
     &           aorespfilename,! in: filename
 
455
     &           npol,          ! in: nr polarization
 
456
     &           nocc,          ! in: nr occupied MOs
 
457
     &           nvirt,         ! in: nr virtual  MOs
 
458
     &           1,             ! in: nr. components
 
459
     &           g_rhs,         ! in: (ntot,3)       GA matrix
 
460
     &           g_rhs_im,      ! in: dummy
 
461
     &           .false.)       ! in: =T if (RE,IM) =F if RE
 
462
 
 
463
c 000000000000 move 2nd subspace to 1st 00000 START
 
464
          shift=ndir
 
465
          m1=shift+1
 
466
          m2=shift+ndir
 
467
          call ga_copy_patch('n',g_rhs,1,ntot,m1,m2, 
 
468
     $                           g_rhs,1,ntot,1 ,ndir)     
 
469
c 000000000000 move 2nd subspace to 1st 00000 END
 
470
c      if (ga_nodeid().eq.0)
 
471
c     & write(*,*) 'FA-BEF get_vecB1_opt2'
 
472
 
 
473
      call get_vecB1_opt2(
 
474
     &                  g_vecB1,    ! out:
 
475
     &                  g_rhs,      ! in : g_rhs vector (occ-virt of g_u)
 
476
     &                  g_u,        ! in : occ-occ of g_u
 
477
     &                  vectors,    ! in : MO vectors
 
478
     &                  nbf,        ! in : nr. basis functions
 
479
     &                  nmo,        ! in : nr. MOs
 
480
     &                  npol,       ! in : nr polarizations
 
481
     &                  nocc,       ! in : nr. occupied MOs
 
482
     &                  nvirt,      ! in : nr. virtual  MOs
 
483
     &                  debug)      ! in : = .true. allow debugging
 
484
 
 
485
c      if (ga_nodeid().eq.0)
 
486
c     & write(*,*) 'FA-AFT get_vecB1_opt2'
 
487
 
 
488
      if (debug) then
 
489
       if (ga_nodeid().eq.0)
 
490
     &  write(*,*) '------- g_vecB1-gb1-nw ---- START'
 
491
       do ispin=1,npol
 
492
        call ga_print(g_vecB1(ispin))
 
493
       enddo
 
494
       if (ga_nodeid().eq.0)
 
495
     &  write(*,*) '------- g_vecB1-gb1-nw ---- END'
 
496
      endif ! end-if-debug
 
497
 
 
498
      do ispin=1,npol
 
499
       if (.not.ga_destroy(g_u(ispin))) call 
 
500
     &   errquit('giao_b1: ga_destroy failed vectors',0,GA_ERR)
 
501
       if (.not.ga_destroy(vectors(ispin))) call 
 
502
     &   errquit('giao_b1: ga_destroy failed vectors',0,GA_ERR)
 
503
      enddo
 
504
      do ispin=1,ndens
 
505
       if (.not.ga_destroy(g_dens(ispin))) call 
597
506
     &    errquit('giao_b1: ga_destroy failed g_dens',0,GA_ERR)
598
 
 
599
 
      if (debug) write (luout,*) 'Perturbed C and Pmat done'
600
 
 
601
 
c     Alll done.
602
 
c     At this point, we don't need to terminate the integrals.
603
 
c     They were terminated by the cphf.
604
 
 
 
507
      enddo ! end-loop-ispin
 
508
c      RHS arrays are no longer needed
 
509
        if (.not.ga_destroy(g_rhs)) call 
 
510
     &     errquit('fiao_f1: ga_destroy failed g_rhs',0,GA_ERR)
605
511
      call ga_sync()
606
 
 
607
512
      return
608
513
 
609
514
 7000 format(/,10x,'B-field perturbed MOs cannot be calculated for',
614
519
     1 /,10x,54(1h-),/,
615
520
     2 10x,'Calculating magnetic field perturbed MO vectors (GIAO)',/,
616
521
     3 10x,54(1h-),/)
617
 
 
618
522
      end
619