← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/property/aor_r1_tensor.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/aor_r1_tensor.F
1
 
 
2
 
      subroutine aor_r1_tensor (rtdb, basis, geom, omega, lstatic, ncomp
3
 
     &   ,  g_smat0, g_dipmag,  g_dipel, g_quadel, g_vectors,
4
 
     &   froct,nbf, nmo, nocct, nvirt, lgiao, lquad, lvelocity,lifetime,
5
 
     &   lmagpert, g_vecE1, g_vecE1_im, alfare, alfaim, betare, betaim)
6
 
 
7
 
c $Id: aor_r1_tensor.F 21176 2011-10-10 06:35:49Z d3y133 $
8
 
      
9
 
c     =================================================================
10
 
      
 
1
      subroutine aor_r1_tensor(
 
2
     &      rtdb,basis,geom,    ! in : handles
 
3
     &      omega,              ! in :
 
4
     &      lstatic,            ! in :
 
5
     &      ncomp,              ! in :  
 
6
     &      g_smat0,            ! in :
 
7
     &      g_dipmag,           ! in : magn -dipole mom AO
 
8
     &      g_dipel,            ! in : elect-dipole mom AO
 
9
     &      g_quadel,           ! in : quadrupole   AO
 
10
     &      g_vectors,          ! in : MOs
 
11
     &      froct,              ! in : set of occupations
 
12
     &      nbf, nmo,           ! in : nr basis, nr MOs
 
13
     &      npol,               ! in : nr. polarizations
 
14
     &      nocct, nvirt,       ! in : nocc,nvirt
 
15
     &      lgiao, lquad,       ! in : logical vars
 
16
     &      lvelocity,          ! in : logical vars
 
17
     &      lifetime,           ! in : logical vars
 
18
     &      lmagpert,           ! in : logical vars
 
19
     &      g_vecE1,g_vecE1_im, ! in : 
 
20
     &      alfare,alfaim,      ! out: electric-electric response matrices
 
21
     &      betare,betaim)      ! out: electric-magnetic response matrices       
 
22
c $Id: aor_r1_tensor.F 23263 2012-12-09 18:38:17Z niri $     
 
23
c     =================================================================   
11
24
c     purpose: calculate linear response tensors
12
 
 
13
25
c     We assume that perturbed MO coefficients have already
14
26
c     been computed elsewhere. 
15
 
 
16
27
c     called from: aoresponse_driver_new
17
 
 
18
28
c     output: alfare, alfaim - field-electric response matrices
19
29
c             betare, betaim - field-magnetic response matrices
20
 
 
21
30
c     =================================================================
 
31
c
 
32
c  Written by J. Autschbach, SUNY Buffalo
 
33
c  Extension to spin-unrestricted case 
 
34
c          by F. Aquino,     Northwestern University 
 
35
c          03-15-12
 
36
c --> Experimental (not published yet)
22
37
 
23
38
      implicit none
24
 
 
25
39
#include "errquit.fh"
26
40
#include "global.fh"
27
41
#include "mafdecls.fh"
33
47
#include "apiP.fh"
34
48
#include "prop.fh"
35
49
#include "bgj.fh"
36
 
 
37
50
c     ---------------------
38
51
c     subroutine arguments:
39
52
c     ---------------------
40
 
 
41
53
      integer rtdb    ! [input] run-time database handle
42
54
      integer basis   ! [input] basis handle
43
55
      integer geom    ! [input] geometry handle
44
 
 
 
56
      integer npol,nocct(npol), nvirt(npol)
 
57
      double precision froct(nbf,npol)
45
58
c     These are all input, too
46
 
      integer g_smat0, g_vectors(2), g_dipel, 
47
 
     &   g_quadel, g_dipmag, g_vecE1(2), g_vecE1_im(2)
 
59
      integer g_smat0, g_vectors(npol), g_dipel, 
 
60
     &        g_quadel, g_dipmag, 
 
61
     &        g_vecE1(2,2),g_vecE1_im(2,2)
48
62
      integer nbf, nmo, ncomp
49
 
      integer nocct(2), nvirt(2)
50
 
      double precision froct(nbf,2)
51
63
      double precision gamwidth, omega
52
64
      logical lgiao, lquad, lvelocity, lifetime, lmagpert, lstatic
53
 
 
 
65
      double precision sum_my ! Added by FA
54
66
c     output:
55
67
      double precision alfare(3,3), alfaim(3,3)
56
68
      double precision betare(3,3), betaim(3,3)
57
 
 
58
 
 
59
69
c     ----------------
60
70
c     local variables:
61
71
c     ----------------
62
 
 
63
 
c     global array handles:
64
 
      
65
 
      integer
66
 
     &   g_work,
67
 
     &   g_temp
68
 
 
 
72
c     global array handles: 
 
73
      integer g_work,g_temp
69
74
c     other local variables: 
70
 
 
71
75
      integer nmot(2), nocvir(2)
72
 
 
73
76
      integer dims(3), chunk(3)
74
77
      integer alo(3), ahi(3), blo(3), bhi(3), clo(3), chi(3)
75
 
 
76
78
      character*(256) cstemp
77
 
 
78
79
      character*(1) direction(3)
79
80
      data direction/'x','y','z'/
80
 
      
81
 
      integer ispin, nspin
82
 
      integer ipm, nocc, nvir, nocv, imo, jmo, nmo1, iresp, idir
 
81
      integer ispin
 
82
      integer ipm, nocc, nvir, nocv, imo, jmo, nmo1, 
 
83
     &        iresp, idir
83
84
      logical debug, dbgmat, 
84
 
     &   lzora, lantisym
 
85
     &        lzora, lantisym
85
86
      double precision sum, scaling
86
87
      double precision tenm8, one, two, zero, half
87
88
      parameter (tenm8=1d-8, one=1d0, two=2d0, zero=0d0, half=one/two)
88
 
 
89
89
c     external functions:
90
 
 
91
 
      double precision ga_trace_diag
92
 
      external ga_trace_diag
93
 
 
 
90
      double precision ga_trace_diag,coeffre,coeffim
 
91
      external ga_trace_diag 
 
92
c ------- Added for unrestricted calc ----- START
 
93
      integer ndir
 
94
      external get_alfaorbeta_reim
 
95
c ------- Added for unrestricted calc ----- END
 
96
      ndir=3 ! nr directions (x,y,z)
94
97
c  ====================================================================
95
 
 
96
 
      debug = .false. .and. ga_nodeid().eq.0 ! .true. during development
 
98
      debug  = .false. .and. ga_nodeid().eq.0 ! .true. during development
97
99
      dbgmat = .false. .and. ga_nodeid().eq.0 ! debug large matrices
98
100
 
 
101
c      debug=.true.
 
102
 
99
103
      if (debug) write (luout,*) 'hello from aor_r1_beta'
100
 
 
101
104
c     the main results are collected in alfare/in, betare/im.
102
105
c     initialize with zeros:
103
 
 
104
 
      alfare = 0
105
 
      alfaim = 0
106
 
      betare = 0
107
 
      betaim = 0
108
 
    
109
 
 
 
106
      do idir=1,ndir
 
107
       do iresp=1,3
 
108
        alfare(idir,iresp) = 0.0d0
 
109
        alfaim(idir,iresp) = 0.0d0
 
110
        betare(idir,iresp) = 0.0d0
 
111
        betaim(idir,iresp) = 0.0d0
 
112
       enddo ! end-loop-iresp
 
113
      enddo ! end-loop-idir
110
114
c     set parameters that control the various computational options
111
115
c     (later we will set most of this by input)
112
 
      nspin      =  1           ! assume closed shell
113
 
      lzora      = .false.      ! not yet available here 
114
 
 
115
 
 
 
116
      lzora  = .false.      ! not yet available here 
116
117
      if (debug) write (luout,*) 'giao, velocity, magpert',
117
118
     &    lgiao, lvelocity, lmagpert
118
 
 
119
119
      lantisym = (lvelocity .or. lmagpert) ! antisymm. perturbation
120
 
 
121
 
 
122
120
c     -----------------------------------------
123
121
c     determine number of occ * virt orbitals
124
122
c     and nmot(1:2) and fix froct, if necessary
125
123
c     -----------------------------------------
126
 
 
127
 
      do ispin = 1,nspin
 
124
      do ispin = 1,npol
128
125
        nocvir(ispin) = nocct(ispin) * nvirt(ispin)
129
126
        nmot(ispin) = nmo
130
127
        if (nmo .lt.nbf) then
132
129
            froct(imo,ispin) = 0d0
133
130
          enddo
134
131
        endif
135
 
      end do
136
 
 
137
 
c     ----------------------------------------------------
138
 
c     start loop over spins (nspin=2 not yet supported)
139
 
c     ----------------------------------------------------
140
 
      
141
 
      do ispin = 1, nspin
142
 
        
143
 
        nmo1 = nmot(ispin)      ! total no.of MOs for this spin
144
 
        nocc = nocct(ispin)     ! occupied MOs
145
 
        nvir = nvirt(ispin)     ! virtual MOs
146
 
        nocv = nocvir(ispin)    ! nocc * nvir
147
 
        
148
 
        
149
 
        
 
132
      enddo ! end-loop-ispin
150
133
c       ------------------------------
151
134
c       allocate some temp. work space
152
135
c       ------------------------------
153
 
        
154
136
        chunk(1) = nbf
155
137
        chunk(2) = -1
156
 
        dims(1) = nbf
157
 
        dims(2) = nbf        
 
138
        dims(1)  = nbf
 
139
        dims(2)  = nbf        
158
140
        write(cstemp,'(a)') 'work'
159
141
        if (.not.nga_create(MT_DBL,2,dims,cstemp(1:4),chunk,
160
142
     &     g_work)) call 
161
143
     &     errquit('aoresponse: nga_create failed: '//cstemp(1:4),
162
144
     &     0,GA_ERR)     
163
145
        call ga_zero (g_work)
164
 
        
 
146
c     ----------------------
 
147
c     start loop over spins
 
148
c     ----------------------
 
149
      if      (npol.eq.1) then
 
150
       coeffre=-2.0d0
 
151
       coeffim=+2.0d0   
 
152
      else if (npol.eq.2) then
 
153
       coeffre=-0.5d0
 
154
       coeffim=+0.5d0  
 
155
      endif
 
156
      do ispin = 1, npol
 
157
        nmo1 = nmot(ispin)      ! total no.of MOs for this spin
 
158
        nocc = nocct(ispin)     ! occupied MOs
 
159
        nvir = nvirt(ispin)     ! virtual MOs
 
160
        nocv = nocvir(ispin)    ! nocc * nvir
165
161
c       allocate intermediate vector for matrix multiplications
166
162
c       used to create the final results
167
 
        
168
163
        write (cstemp,'(a)') 'aor_beta: temp1'
169
164
        if(.not.ga_create(MT_DBL, nbf, nocc, trim(cstemp),
170
165
     &     -1,-1,g_temp))
171
166
     &     call errquit (trim(cstemp),0,GA_ERR)
172
167
        if (debug) write (luout,*) 'g_temp allocated'
173
 
 
174
168
c       -------------------------------------------------------
175
169
c       (A) calculate optical rotation beta from C(E) H(B) C(0)
176
170
c       ------------------------------------------------------
177
 
        
178
 
        
179
 
        
180
171
c       ---------------------------------------------------------
181
172
c       solution of CPKS is in g_vecE1. For the OR and a length-gauge
182
173
c       dipole response (E-field perturbation) we need the difference of
185
176
c       perturbing field has an antisymmetric AO matrix and 
186
177
c       we have to ADD the vectors instead. 
187
178
c       ---------------------------------------------------------
188
 
        
 
179
        if (debug) then
 
180
         if (ga_nodeid().eq.0) then
 
181
          write(*,100) ncomp,lantisym,lifetime
 
182
 100      format('(comp,lantisym,lifetime)=(',i3,',',L1,',',L1,')')
 
183
         endif
 
184
         if (ga_nodeid().eq.0)
 
185
     &    write(*,*) '--g_vecE1-BEF-IF ------ START'
 
186
          call ga_print(g_vecE1(ispin,1))
 
187
         if (ga_nodeid().eq.0)
 
188
     &    write(*,*) '--g_vecE1-BEF-IF ------ END'
 
189
         if (lifetime) then
 
190
          if (ga_nodeid().eq.0)
 
191
     &      write(*,*) '--g_vecE1-im-BEF-IF ------ START'
 
192
            call ga_print(g_vecE1_im(ispin,1))
 
193
          if (ga_nodeid().eq.0)
 
194
     &      write(*,*) '--g_vecE1-im-BEF-IF ------ END'
 
195
         endif ! end-if-lifetime
 
196
        endif ! end-if-debug
189
197
        if (ncomp.gt.1) then
190
198
          if (lantisym) then
191
 
            call ga_add(1d0, g_vecE1(1), 1d0,  g_vecE1(2),
192
 
     &         g_vecE1(1))
 
199
            call ga_add(1d0, g_vecE1(ispin,1), 
 
200
     &                  1d0, g_vecE1(ispin,2),
 
201
     &                       g_vecE1(ispin,1))
193
202
            if (lifetime) then
194
 
              call ga_add(1d0, g_vecE1_im(1), 1d0,  g_vecE1_im(2),
195
 
     &           g_vecE1_im(1))
 
203
              call ga_add(1d0, g_vecE1_im(ispin,1), 
 
204
     &                    1d0, g_vecE1_im(ispin,2),
 
205
     &                         g_vecE1_im(ispin,1))
196
206
            end if              ! lifetime
197
207
          else                  ! lantisym ? 
198
 
            call ga_add(1d0, g_vecE1(1), -1d0,  g_vecE1(2),
199
 
     &         g_vecE1(1))
 
208
            call ga_add(1d0, g_vecE1(ispin,1), 
 
209
     &                 -1d0, g_vecE1(ispin,2),
 
210
     &                       g_vecE1(ispin,1))
200
211
            if (lifetime) then
201
 
              call ga_add(1d0, g_vecE1_im(1), -1d0,  g_vecE1_im(2),
202
 
     &           g_vecE1_im(1))
 
212
              call ga_add(1d0, g_vecE1_im(ispin,1),
 
213
     &                   -1d0, g_vecE1_im(ispin,2),
 
214
     &                         g_vecE1_im(ispin,1))
203
215
            end if              ! lifetime
204
216
          end if                ! lantisym
205
217
        endif                   ! ncomp.gt.1
206
 
               
207
 
        do idir = 1,3           ! direction of the perturbing field
208
 
                    
209
 
          do iresp = 1,3
210
 
                                                 
211
 
            alo(1) = 1
212
 
            ahi(1) = nbf
213
 
            alo(2) = 1
214
 
            ahi(2) = nbf
215
 
            alo(3) = iresp      ! pick the correct B-field direction
216
 
            ahi(3) = iresp
217
 
            blo(1) = 1
218
 
            bhi(1) = nbf
219
 
            blo(2) = 1
220
 
            bhi(2) = nocc
221
 
            clo(1) = 1
222
 
            chi(1) = nbf
223
 
            clo(2) = 1
224
 
            chi(2) = nocc
225
 
                        
226
 
            call ga_zero(g_temp)
227
 
            call nga_matmul_patch('n','n',1d0,0d0,
228
 
     &         g_dipmag,alo,ahi,
229
 
     &         g_vectors(ispin),blo,bhi,
230
 
     &         g_temp,clo,chi)
231
 
            
232
 
            if (debug) write (luout,*)
233
 
     &         'beta: H(B) C(0) intermediate complete'
234
 
            
235
 
            alo(1) = 1
236
 
            ahi(1) = nocc
237
 
            alo(2) = 1
238
 
            ahi(2) = nbf
239
 
            alo(3) = idir
240
 
            ahi(3) = idir
241
 
            blo(1) = 1
242
 
            bhi(1) = nbf
243
 
            blo(2) = 1
244
 
            bhi(2) = nocc
245
 
            clo(1) = 1
246
 
            chi(1) = nocc
247
 
            clo(2) = 1
248
 
            chi(2) = nocc
249
 
            
250
 
            call ga_zero(g_work)
251
 
            call nga_matmul_patch('t','n',1d0,0d0,
252
 
     &         g_vecE1,alo,ahi,
253
 
     &         g_temp,blo,bhi,
254
 
     &         g_work,clo,chi)
255
 
                        
256
 
            
257
 
c           the factor of two is for the orbital occupations,
258
 
c           assuming that ispin is never equal to two
259
 
            
260
 
            sum = 2d0 * ga_trace_diag(g_work)
261
 
            
262
 
c           the minus sign is consistent with the original aoresponse
263
 
c           routine. We have to use - to get consistent results
264
 
c           with the output from aor_r1_beta.F
265
 
            betare(idir,iresp) = betare(idir,iresp) - sum
266
 
 
267
 
 
268
 
            if (lifetime) then
269
 
              call ga_zero(g_work)
270
 
              call nga_matmul_patch('t','n',1d0,0d0,
271
 
     &           g_vecE1_im,alo,ahi,
272
 
     &           g_temp,blo,bhi,
273
 
     &           g_work,clo,chi)
274
 
              
275
 
            
276
 
            sum = 2d0 * ga_trace_diag(g_work)
277
 
            betaim(idir,iresp) = betaim(idir,iresp) + sum
278
 
 
279
 
          end if                ! lifetime
280
 
          
281
 
          if (debug) write (luout,*)
282
 
     &       'beta: C(E) H(B) C(0) complete'
283
 
 
284
 
          end do                ! iresp
285
 
        end do                  ! idir
286
 
            
287
 
            
 
218
        if (debug) then
 
219
         if (ga_nodeid().eq.0)
 
220
     &    write(*,*) '--g_vecE1-AFT-IF ------ START'
 
221
          call ga_print(g_vecE1(ispin,1))
 
222
         if (ga_nodeid().eq.0)
 
223
     &    write(*,*) '--g_vecE1-AFT-IF ------ END'
 
224
         if (lifetime) then
 
225
          if (ga_nodeid().eq.0)
 
226
     &     write(*,*) '--g_vecE1-im-AFT-IF ------ START'
 
227
           call ga_print(g_vecE1_im(ispin,1))
 
228
          if (ga_nodeid().eq.0)
 
229
     &     write(*,*) '--g_vecE1-im-AFT-IF ------ END'
 
230
         endif ! end-if-lifetime
 
231
        endif ! end-if-debug
 
232
c       if (ga_nodeid().eq.0)
 
233
c     &  write(*,*) '==== compute beta ====== START'
 
234
       ndir=3 ! nr. directions (x,y,z)
 
235
       do idir = 1,ndir        ! direction of the perturbing field
 
236
        do iresp = 1,3
 
237
 
 
238
         call get_alfaorbeta_reim(
 
239
     &            betare(idir,iresp), ! in/out: alpha or beta real part
 
240
     &            betaim(idir,iresp), ! in/out: alpha or beta im   part
 
241
     &            g_vecE1(ispin,1),   ! in : 1st-order pert vec RE
 
242
     &            g_vecE1_im(ispin,1),! in : 1st-order pert vec IM
 
243
     &            g_dipmag,           ! in : dipole electric or magnetic
 
244
     &            g_vectors(ispin),   ! in : MO vectors
 
245
     &            idir,               ! in : = 1,2,3=x,y,z directions
 
246
     &            iresp,              ! in : = 1,2,3
 
247
     &            coeffre,coeffim,1,  ! in : (coeffre,coeffim,caseAO)
 
248
     &            nbf,                ! in : nr. basis functions
 
249
     &            nocc,               ! in : nr. occupied alpha or beta
 
250
     &            lifetime,           ! in : logical var for damping
 
251
     &            debug,              ! in : logical var for debugging
 
252
     &            g_temp)             ! in : scratch GA array
 
253
 
 
254
         if (debug) then 
 
255
          if (ga_nodeid().eq.0) then
 
256
           write(*,1) ispin,idir,iresp,
 
257
     &                betare(idir,iresp),betaim(idir,iresp)
 
258
 1         format('beta(',i3,',',i3,',',i3,
 
259
     &            ')=(',f15.8,',',f15.8,')')
 
260
          endif
 
261
         endif ! end-if-debug
 
262
        enddo  ! end-loop-iresp (responding field components)
 
263
       enddo ! end-loop-idir    (perturbing E-field components)
 
264
c       if (ga_nodeid().eq.0)
 
265
c     &  write(*,*) '==== compute beta ====== END'
288
266
c       --------------------------------------
289
267
c       (B) calculate alfa from C(E) h(E) C(0)
290
268
c       --------------------------------------
291
 
 
292
269
c       --------------------------------------------------------- 
293
270
c       For alfa we need the sum of the +/- components no matter
294
271
c       if we use
295
272
c       length or velocity gauge so we add twice the icomp=2 component
296
273
c       back into g_vecE1(1)
297
274
c       ---------------------------------------------------------
298
 
        
299
275
        if (ncomp.gt.1) then
300
276
          if (lantisym) then
301
277
            continue
302
278
          else
303
 
            call ga_add(1d0, g_vecE1(1), 2d0,  g_vecE1(2),
304
 
     &         g_vecE1(1))
 
279
            call ga_add(1d0, g_vecE1(ispin,1),
 
280
     &                  2d0, g_vecE1(ispin,2),
 
281
     &                       g_vecE1(ispin,1))
305
282
            if (lifetime) then
306
 
              call ga_add(1d0, g_vecE1_im(1), 2d0,  g_vecE1_im(2),
307
 
     &           g_vecE1_im(1))
 
283
              call ga_add(1d0, g_vecE1_im(ispin,1), 
 
284
     &                    2d0, g_vecE1_im(ispin,2),
 
285
     &                         g_vecE1_im(ispin,1))
308
286
            end if              ! lifetime
309
287
          end if                ! lantisym
310
288
        endif                   ! ncomp.gt.1
311
 
 
312
 
        do idir = 1,3
313
 
          do iresp = 1,3
314
 
          
315
 
            alo(1) = 1
316
 
            ahi(1) = nbf
317
 
            alo(2) = 1
318
 
            ahi(2) = nbf
319
 
            alo(3) = iresp      ! pick direction iresp for g_dipel
320
 
            ahi(3) = iresp
321
 
            blo(1) = 1
322
 
            bhi(1) = nbf
323
 
            blo(2) = 1
324
 
            bhi(2) = nocc
325
 
            clo(1) = 1
326
 
            chi(1) = nbf
327
 
            clo(2) = 1
328
 
            chi(2) = nocc
329
 
            
330
 
            
331
 
            call ga_zero(g_temp)
332
 
            call nga_matmul_patch('n','n',1d0,0d0,
333
 
     &         g_dipel,alo,ahi,
334
 
     &         g_vectors(ispin),blo,bhi,
335
 
     &         g_temp,clo,chi)
336
 
            
337
 
            if (debug) write (luout,*)
338
 
     &         'alfa: h(E) C(0) intermediate complete'
339
 
            
340
 
            
341
 
            alo(1) = 1
342
 
            ahi(1) = nocc
343
 
            alo(2) = 1
344
 
            ahi(2) = nbf
345
 
            alo(3) = idir
346
 
            ahi(3) = idir
347
 
            blo(1) = 1
348
 
            bhi(1) = nbf
349
 
            blo(2) = 1
350
 
            bhi(2) = nocc
351
 
            clo(1) = 1
352
 
            chi(1) = nocc
353
 
            clo(2) = 1
354
 
            chi(2) = nocc
355
 
            
356
 
            call ga_zero(g_work)
357
 
            
358
 
            call nga_matmul_patch('t','n',1d0,0d0,
359
 
     &         g_vecE1,alo,ahi,
360
 
     &         g_temp,blo,bhi,
361
 
     &         g_work,clo,chi)
362
 
            
363
 
                       
364
 
c           the factor of two is for the orbital occupations,
365
 
c           assuming that ispin is never equal to two
366
 
            
367
 
            sum = 2d0 * ga_trace_diag(g_work)
368
 
            
369
 
            alfare(idir,iresp) = alfare(idir,iresp) - sum
370
 
 
371
 
            if (lifetime) then
372
 
 
373
 
              call ga_zero(g_work)
374
 
              
375
 
              call nga_matmul_patch('t','n',1d0,0d0,
376
 
     &           g_vecE1_im,alo,ahi,
377
 
     &           g_temp,blo,bhi,
378
 
     &           g_work,clo,chi)
379
 
              
380
 
              
381
 
              sum = 2d0 * ga_trace_diag(g_work)
382
 
              
383
 
              alfaim(idir,iresp) =  alfaim(idir,iresp) + sum
384
 
            end if              ! lifetime
385
 
            
386
 
            if (debug) write (luout,*) 'alfa C(E) h(E) C(0) complete'
387
 
                        
388
 
          enddo                 ! iresp = 1,3       
389
 
 
390
 
c         -----------------------------------------
391
 
c         end loop over responding field components
392
 
c         -----------------------------------------
393
 
 
394
 
          
395
 
        end do                  ! idir = 1,3
396
 
 
397
 
c       -------------------------------------------
398
 
c       end loop over perturbing E-field components
399
 
c       -------------------------------------------
400
 
 
401
 
c       -----------------
402
 
c       deallocate memory
403
 
c       -----------------
404
 
 
 
289
c       if (ga_nodeid().eq.0)
 
290
c     &  write(*,*) '==== compute alfa ====== START'
 
291
       ndir=3 ! nr. directions (x,y,z)
 
292
       do idir = 1,ndir        ! direction of the perturbing field
 
293
        do iresp = 1,3
 
294
 
 
295
         call get_alfaorbeta_reim(
 
296
     &            alfare(idir,iresp), ! out: alpha or beta real part
 
297
     &            alfaim(idir,iresp), ! out: alpha or beta im   part
 
298
     &            g_vecE1(ispin,1),   ! in : 1st-order pert vec RE
 
299
     &            g_vecE1_im(ispin,1),! in : 1st-order pert vec IM
 
300
     &            g_dipel,            ! in : dipole electric or magnetic
 
301
     &            g_vectors(ispin),   ! in : MO vectors
 
302
     &            idir,               ! in : = 1,2,3=x,y,z directions
 
303
     &            iresp,              ! in : = 1,2,3
 
304
     &            coeffre,coeffim,1,  ! in : (coeffre,coeffim,caseAO)
 
305
     &            nbf,                ! in : nr. basis functions
 
306
     &            nocc,               ! in : nr. occupied alpha or beta
 
307
     &            lifetime,           ! in : logical var for damping
 
308
     &            debug,              ! in : logical var for debugging
 
309
     &            g_temp)             ! in : scratch GA array
 
310
 
 
311
         if (debug) then 
 
312
          if (ga_nodeid().eq.0) then
 
313
           write(*,2) ispin,idir,iresp,
 
314
     &                alfare(idir,iresp),alfaim(idir,iresp)
 
315
 2         format('alfa(',i3,',',i3,',',i3,
 
316
     &            ')=(',f15.8,',',f15.8,')')
 
317
          endif
 
318
         endif ! end-if-debug
 
319
        enddo  ! end-loop-iresp (responding field components)
 
320
       enddo ! end-loop-idir    (perturbing E-field components)
 
321
c       if (ga_nodeid().eq.0)
 
322
c     &  write(*,*) '==== compute alfa ====== END'
 
323
c ============ visualize-1 (alfa,beta) ========== START
 
324
c        if (ga_nodeid().eq.0) then
 
325
c         do idir=1,ndir
 
326
c          do iresp=1,3
 
327
c           write(*,10) ispin,idir,iresp,
 
328
c     &       alfare(idir,iresp),alfaim(idir,iresp),
 
329
c     &       betare(idir,iresp),betaim(idir,iresp)
 
330
c 10       format('FA-tensor:(ispin,idir,iresp)=(',
 
331
c     &            i3,',',i3',',i3,')',
 
332
c     &           ' alfa(re,im)=(',f15.8,',',f15.8,')',
 
333
c     &           ' beta(re,im)=(',f15.8,',',f15.8,')')
 
334
c          enddo ! end-loop-iresp
 
335
c         enddo ! end-loop-idir
 
336
c        endif
 
337
c ============ visualize-1 (alfa,beta) ========== END
405
338
          if (.not.ga_destroy(g_temp))
406
339
     &       call errquit
407
340
     &       ('aor_beta: ga_destroy failed g_temp',
408
341
     &       0,GA_ERR)
409
 
 
 
342
      enddo ! end-loop-ispin
 
343
c       -----------------
 
344
c       deallocate memory
 
345
c       -----------------
410
346
 
411
347
        if (.not.ga_destroy(g_work))
412
348
     &     call 
413
349
     &     errquit('aoresponse: ga_destroy failed g_work',
414
350
     &     0,GA_ERR)
415
 
 
416
 
        
417
 
      enddo                     ! ispin = 1,2 from way above
418
 
              
419
 
c     ---------------------------------------------------------------
420
 
c     end loop over spin components (which we don't use right now
421
 
c     since nspin is forced to be 1 at the beginning of this routine)
422
 
c     ---------------------------------------------------------------
423
 
            
424
 
 
425
351
c     it seems that if we use GIAOs everything is off by a factor of
426
352
c     two, so we need to scale betare, betaim. If we have static
427
353
c     response then there is a factor of two missing everywhere
428
354
c     because we don't add C(+) and C(-) for the electric field.
429
355
 
 
356
c       if (ga_nodeid().eq.0)
 
357
c     &  write(*,*) 'FA-lgiao=',lgiao
 
358
 
430
359
      if (lgiao) then
431
360
        scaling = half
432
 
        do idir = 1,3
 
361
        do idir = 1,ndir ! direction of the perturbing field (x,y,z)
433
362
          do iresp = 1,3
434
363
            betare(idir, iresp) = betare(idir, iresp) * scaling
435
364
            betaim(idir, iresp) = betaim(idir, iresp) * scaling
436
365
          end do
437
366
        end do
438
 
      end if                    ! lstatic
 
367
      end if                    ! lgiao
 
368
c ============ visualize-1 (alfa,beta) ========== START
 
369
c        if (ga_nodeid().eq.0) then
 
370
c         do idir=1,ndir
 
371
c          do iresp=1,3
 
372
c           write(*,11) ispin,idir,iresp,
 
373
c     &       alfare(idir,iresp),alfaim(idir,iresp),
 
374
c     &       betare(idir,iresp),betaim(idir,iresp)
 
375
c 11       format('FA-AFT-lgiao:(ispin,idir,iresp)=(',
 
376
c     &            i3,',',i3',',i3,')',
 
377
c     &           ' alfa(re,im)=(',f15.8,',',f15.8,')',
 
378
c     &           ' beta(re,im)=(',f15.8,',',f15.8,')')
 
379
c          enddo ! end-loop-iresp
 
380
c         enddo ! end-loop-idir
 
381
c        endif
 
382
c ============ visualize-1 (alfa,beta) ========== END
 
383
c       if (ga_nodeid().eq.0)
 
384
c     &  write(*,*) 'FA-lstatic=',lstatic
439
385
 
440
386
      if (lstatic) then
441
387
        scaling = two
442
 
        do idir = 1,3
 
388
        do idir = 1,ndir ! direction of the perturbing field (x,y,z)
443
389
          do iresp = 1,3
444
390
            alfare(idir, iresp) = alfare(idir, iresp) * scaling
445
391
            alfaim(idir, iresp) = alfaim(idir, iresp) * scaling
449
395
        end do
450
396
      end if                    ! lstatic
451
397
 
 
398
c ============ visualize-1 (alfa,beta) ========== START
 
399
c        if (ga_nodeid().eq.0) then
 
400
c         do idir=1,ndir
 
401
c          do iresp=1,3
 
402
c           write(*,12) ispin,idir,iresp,
 
403
c     &       alfare(idir,iresp),alfaim(idir,iresp),
 
404
c     &       betare(idir,iresp),betaim(idir,iresp)
 
405
c 12       format('FA-AFT-lstatic:(ispin,idir,iresp)=(',
 
406
c     &            i3,',',i3',',i3,')',
 
407
c     &           ' alfa(re,im)=(',f15.8,',',f15.8,')',
 
408
c     &           ' beta(re,im)=(',f15.8,',',f15.8,')')
 
409
c          enddo ! end-loop-iresp
 
410
c         enddo ! end-loop-idir
 
411
c        endif
 
412
c ============ visualize-1 (alfa,beta) ========== END
452
413
 
453
414
c     for magnetic perturbation, alfa is the G' tensor, and
454
415
c     beta is the paramagnetic mag-mag response. 
457
418
c     omega we want the beta tensor, not G'. 
458
419
c     also consider the case of velocity formalisms where 
459
420
c     additional divisions by omega occur. 
460
 
        
461
 
      do idir = 1,3
 
421
 
 
422
c      if (ga_nodeid().eq.0) then
 
423
c       write(*,14) lmagpert,lstatic,lvelocity,omega
 
424
c 14    format('(lmagpert,lstatic,lvelocity,omega)=(',
 
425
c     &        L1,',',L1,',',L1,',',f15.8,')')
 
426
c      endif
 
427
 
 
428
      do idir = 1,ndir ! direction of the perturbing field (x,y,z)
462
429
        do iresp = 1,3
463
430
          if (lmagpert .and. .not.lstatic) then
464
431
c           case I: magnetic perturbation. alpha = G', beta = Chi-p
465
 
            alfare(idir,iresp) = alfare(idir,iresp) / omega
466
 
            alfaim(idir,iresp) = alfaim(idir,iresp) / omega
 
432
c            if (ga_nodeid().eq.0)
 
433
c     &       write(*,*) 'FA-enter-caseI'
 
434
              alfare(idir,iresp) = alfare(idir,iresp) / omega
 
435
              alfaim(idir,iresp) = alfaim(idir,iresp) / omega
467
436
            if (lvelocity) then
468
437
              alfare(idir,iresp) = alfare(idir,iresp) / omega
469
438
              alfaim(idir,iresp) = alfaim(idir,iresp) / omega
470
439
            endif
471
440
          else if (.not.lmagpert .and. .not.lstatic) then
 
441
c            if (ga_nodeid().eq.0)
 
442
c     &       write(*,*) 'FA-enter-caseII'
472
443
c           case II: alpha = polarizability, beta = G'
473
 
            betare(idir,iresp) = betare(idir,iresp) / omega
474
 
            betaim(idir,iresp) = betaim(idir,iresp) / omega
 
444
              betare(idir,iresp) = betare(idir,iresp) / omega
 
445
              betaim(idir,iresp) = betaim(idir,iresp) / omega
475
446
            if (lvelocity) then
 
447
c              if (ga_nodeid().eq.0) 
 
448
c     &         write(*,*) 'FA-enter-caseII-lvelocity'
476
449
              betare(idir,iresp) = betare(idir,iresp) / omega
477
450
              alfare(idir,iresp) = alfare(idir,iresp) / omega**2
478
451
              betaim(idir,iresp) = betaim(idir,iresp) / omega
481
454
          endif
482
455
        enddo
483
456
      enddo
484
 
 
485
 
     
 
457
c ============ visualize-1 (alfa,beta) ========== START
 
458
c        if (ga_nodeid().eq.0) then
 
459
c         do idir=1,ndir
 
460
c          do iresp=1,3
 
461
c           write(*,13) ispin,idir,iresp,
 
462
c     &       alfare(idir,iresp),alfaim(idir,iresp),
 
463
c     &       betare(idir,iresp),betaim(idir,iresp)
 
464
c 13       format('FA-AFT-lmagpert:(ispin,idir,iresp)=(',
 
465
c     &            i3,',',i3',',i3,')',
 
466
c     &           ' alfa(re,im)=(',f15.8,',',f15.8,')',
 
467
c     &           ' beta(re,im)=(',f15.8,',',f15.8,')')
 
468
c          enddo ! end-loop-iresp
 
469
c         enddo ! end-loop-idir
 
470
c        endif
 
471
c ============ visualize-1 (alfa,beta) ========== END
486
472
c     ----------------
487
473
c     all done. return
488
474
c     ----------------
489
 
                  
490
 
      
491
475
c     ==================================================================
492
 
      
493
476
      return
494
 
      
495
477
      end
496