← 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_beta_anl_tools.F

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

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/property/aor_r1_beta_anl_tools.F
 
1
      subroutine debug_update_g_anl_quad1(
 
2
     &            g_quadel,
 
3
     &            g_anl, 
 
4
     &            g_vecE1,
 
5
     &            g_vectors,
 
6
     &            ispin)
 
7
c Purpose: Debugging input to update_g_anl_quad
 
8
c
 
9
c Author : Fredy W. Aquino
 
10
c Date   : 03-15-12
 
11
c Note.- Modified from original aoresponse source code
 
12
c        for extension to spin-unrestricted case
 
13
c        original aoresponse source code was written by 
 
14
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
15
c --> Experimental (not published yet)
 
16
        
 
17
      implicit none
 
18
#include "errquit.fh"
 
19
#include "global.fh"
 
20
#include "mafdecls.fh"
 
21
#include "msgids.fh"
 
22
#include "stdio.fh"
 
23
       integer ispin
 
24
       integer g_quadel,g_anl, 
 
25
     &         g_vecE1(2,2),g_vectors(2)
 
26
           if (ga_nodeid().eq.0)
 
27
     &     write(*,*) '---- g_quadel-------- START'
 
28
           call ga_print(g_quadel)
 
29
           if (ga_nodeid().eq.0)
 
30
     &     write(*,*) '---- g_quadel-------- END'
 
31
 
 
32
           if (ga_nodeid().eq.0) then
 
33
            write(*,17) ispin
 
34
 17         format('---- g_anl-lquad-bef(',i3,')-------- START')
 
35
           endif
 
36
           call ga_print(g_anl)
 
37
           if (ga_nodeid().eq.0) then
 
38
            write(*,18) ispin
 
39
 18         format('---- g_anl-lquad-bef(',i3,')-------- END')
 
40
           endif
 
41
 
 
42
           if (ga_nodeid().eq.0) then
 
43
            write(*,13) ispin
 
44
 13         format('---- g_vecE1(',i3,')-------- START')
 
45
           endif
 
46
           call ga_print(g_vecE1(ispin,1))
 
47
           if (ga_nodeid().eq.0) then
 
48
            write(*,14) ispin
 
49
 14         format('---- g_vecE1(',i3,')-------- END')
 
50
           endif
 
51
 
 
52
           if (ga_nodeid().eq.0) then
 
53
            write(*,15) ispin
 
54
 15         format('---- g_vectors(',i3,')-------- START')
 
55
           endif
 
56
           call ga_print(g_vectors(ispin))
 
57
           if (ga_nodeid().eq.0) then
 
58
            write(*,16) ispin
 
59
 16         format('---- g_vectors(',i3,')-------- END')
 
60
           endif
 
61
      return
 
62
      end
 
63
 
 
64
      subroutine update_g_anl_quad(
 
65
     &            g_anl,    ! in/out:
 
66
     &            g_quadel, ! in    : quadrupole AO integral
 
67
     &            g_vecE1,  ! in    : 1st-order elect. pert. MO vector
 
68
     &            g_vectors,! in    : MO vector
 
69
     &            scaling,  ! in    : scaling factor
 
70
     &            g_temp,   ! in    : scratch GA array
 
71
     &            g_work,   ! in    : scratch GA array
 
72
     &            idir,     ! in    : =1,2,3=x,y,z
 
73
     &            nocc,     ! in    : nr. occupied MOs      
 
74
     &            nbf,      ! in    : nr. basis functions  
 
75
     &            debug,    ! in    : debugging flag  
 
76
     &            lstatic)  ! in    : static    flag
 
77
c
 
78
c Author : Fredy W. Aquino
 
79
c Date   : 03-15-12
 
80
c Note.- Modified from original aoresponse source code
 
81
c        for extension to spin-unrestricted case
 
82
c        original aoresponse source code was written by 
 
83
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
84
c --> Experimental (not published yet)
 
85
        
 
86
      implicit none
 
87
#include "errquit.fh"
 
88
#include "global.fh"
 
89
#include "mafdecls.fh"
 
90
#include "msgids.fh"
 
91
#include "geom.fh"
 
92
#include "rtdb.fh"
 
93
#include "bas.fh"
 
94
#include "stdio.fh"
 
95
#include "prop.fh"
 
96
      integer npol
 
97
      integer g_anl,   ! in/out
 
98
     &        g_quadel,  
 
99
     &        g_vecE1,g_vectors,
 
100
     &        g_temp,g_work   
 
101
      logical lstatic,debug
 
102
      integer nbf,nocc,iresp,idir
 
103
      integer alo(3),ahi(3), 
 
104
     &        blo(3),bhi(3), 
 
105
     &        clo(3),chi(3)      
 
106
      integer ind_dkl(2,3)
 
107
      data ind_dkl/ 2, 3,  ! dkl=123
 
108
     &              3, 1,  ! dkl=231
 
109
     &              1, 2 / ! dkl=312
 
110
      integer ind_p(2,2),p,ind1,ind2,k,l
 
111
      data ind_p/1, 2,   ! p12   p=1
 
112
     &           2, 1 /  ! p21   p=2
 
113
      integer LCred(2)
 
114
      data LCred/1,-1/
 
115
      integer qindex(3,3)
 
116
      double precision scl,scaling,one,two,four
 
117
      parameter(one=1.0d0,two=2.0d0,four=4.0d0)
 
118
      external get_C1MC
 
119
c     define translation table for quadrupole indices in
 
120
c     packed storage
 
121
c     XX=1, XY=YX=2, XZ=ZX=3, YY=4, YZ=ZY=5, ZZ=6
 
122
      qindex(1,1) = 1
 
123
      qindex(1,2) = 2
 
124
      qindex(2,1) = 2
 
125
      qindex(1,3) = 3
 
126
      qindex(3,1) = 3
 
127
      qindex(2,2) = 4
 
128
      qindex(2,3) = 5
 
129
      qindex(3,2) = 5
 
130
      qindex(3,3) = 6             
 
131
      do p=1,2 
 
132
        ind1=ind_p(1,p)
 
133
        ind2=ind_p(2,p)
 
134
c Allowing: (dkl)=(123,132,231,213,312,321)  d=idir
 
135
        k=ind_dkl(ind1,idir)  
 
136
        l=ind_dkl(ind2,idir)   
 
137
        scl=-LCred(p)*scaling               
 
138
        iresp = qindex(l,idir)  
 
139
 
 
140
        call get_C1MC(g_work,    ! out: C(E) M C
 
141
     &                g_vecE1,   ! in : 1st-order pert vec RE
 
142
     &                g_quadel,  ! in : dipole electric or magnetic
 
143
     &                g_vectors, ! in : MO vectors
 
144
     &                k,         ! in : = 1,2,3=x,y,z directions
 
145
     &                iresp,     ! in : = 1,2,3
 
146
     &                1,         ! in : indices in (alo,ahi)(3) (blo,bhi)(3)
 
147
     &                nbf,       ! in : nr. basis functions
 
148
     &                nocc,      ! in : nr. occupied alpha or beta
 
149
     &                debug,     ! in : logical var for debugging
 
150
     &                g_temp)    ! in : scratch GA array -> out
 
151
 
 
152
        clo(1) = 1
 
153
        chi(1) = nocc
 
154
        clo(2) = 1
 
155
        chi(2) = nocc           
 
156
        call nga_add_patch(scl,g_work,clo,chi,
 
157
     &                     one,g_anl ,clo,chi,
 
158
     &                         g_anl ,clo,chi)  
 
159
      enddo ! end-loop-p  
 
160
      return
 
161
      end
 
162
 
 
163
      subroutine update_g_anl_elmag(
 
164
     &              g_anl,  ! in/out:
 
165
     &              g_vecE1,! in    : C(E) pert MO vect.
 
166
     &              g_smatx,! in    : 0th/1st overlap deriv.
 
167
     &              g_vec,  ! in    : MO vect or C(B) MO vect
 
168
     &              scaling,! in    : scaling factor
 
169
     &              idir,   ! in    : direction (x,y or z)
 
170
     &              caseAO, ! in    : 3 or 4
 
171
     &              nbf,    ! in    : nr. basis functions
 
172
     &              nocc,   ! in    : nr. occupied MOs
 
173
     &              lstatic,! in    : flag for static calc.
 
174
     &              debug,  ! in    : flag for debugging
 
175
     &              g_work, ! in    : scratch GA array
 
176
     &              g_temp) ! in    : scratch GA array
 
177
c
 
178
c Author : Fredy W. Aquino
 
179
c Date   : 03-15-12
 
180
c Note.- Modified from original aoresponse source code
 
181
c        for extension to spin-unrestricted case
 
182
c        original aoresponse source code was written by 
 
183
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
184
c --> Experimental (not published yet)
 
185
        
 
186
      implicit none
 
187
#include "errquit.fh"
 
188
#include "global.fh"
 
189
#include "mafdecls.fh"
 
190
#include "msgids.fh"
 
191
#include "bas.fh"
 
192
#include "stdio.fh"
 
193
       integer clo(3),chi(3)
 
194
       integer caseAO,nbf,nocc,idir
 
195
       logical lstatic,debug
 
196
       integer g_smatx,g_vecE1,g_vec,g_anl,
 
197
     &         g_temp,g_work ! scratch GA arrays
 
198
       double precision scaling,one,two,four
 
199
       parameter (one=1d0,two=2d0,four=4.0d0)
 
200
       external get_C1MC
 
201
       if (caseAO.ne.3 .and. caseAO.ne.4) then
 
202
        call errquit('update_g_anl_elmag: caseAO ne 3 or 4',
 
203
     &               1,INPUT_ERR)
 
204
       endif
 
205
 
 
206
         call get_C1MC(
 
207
     &           g_work,  ! out: C(E) M C
 
208
     &           g_vecE1, ! in : 1st-order pert vec RE
 
209
     &           g_smatx, ! in : dipole electric or magnetic
 
210
     &           g_vec,   ! in : MO vectors
 
211
     &           idir,    ! in : = 1,2,3=x,y,z directions
 
212
     &           0,       ! in : for caseAO=3,4 is dummy
 
213
     &           caseAO,  ! in : caseAO:(ind1,ind2)=(idir,1)
 
214
     &           nbf,     ! in : nr. basis functions
 
215
     &           nocc,    ! in : nr. occupied alpha or beta
 
216
     &           debug,   ! in : logical var for debugging
 
217
     &           g_temp)  ! in : scratch GA array   
 
218
 
 
219
       if (debug) then 
 
220
        if (ga_nodeid().eq.0) then
 
221
         write(*,5) 1
 
222
 5       format('---- g_work(',i3,')-------- START')
 
223
        endif
 
224
        call ga_print(g_work)
 
225
        if (ga_nodeid().eq.0) then
 
226
         write(*,6) 1
 
227
 6       format('---- g_work(',i3,')-------- END')
 
228
        endif
 
229
       endif ! end-if-debug
 
230
       clo(1) = 1
 
231
       chi(1) = nocc
 
232
       clo(2) = 1
 
233
       chi(2) = nocc
 
234
       call nga_add_patch(scaling,g_work,clo,chi,
 
235
     &                        one,g_anl ,clo,chi,
 
236
     &                            g_anl ,clo,chi)
 
237
      return
 
238
      end
 
239
 
 
240
      subroutine get_trace_ganl(g_anl,  ! in/out:
 
241
     &                          g_tran, ! in:
 
242
     &                          idir,   ! in: 1,2,3=x,y,z
 
243
     &                          oprint, ! in: logical var
 
244
     &                          nocc)   ! in: nr. occupied MOs
 
245
c
 
246
c Author : Fredy W. Aquino
 
247
c Date   : 03-15-12
 
248
c Note.- Modified from original aoresponse source code
 
249
c        for extension to spin-unrestricted case
 
250
c        original aoresponse source code was written by 
 
251
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
252
c --> Experimental (not published yet)
 
253
        
 
254
      implicit none
 
255
#include "errquit.fh"
 
256
#include "global.fh"
 
257
#include "mafdecls.fh"
 
258
#include "msgids.fh"
 
259
#include "stdio.fh"
 
260
      integer nocc,i,idir,
 
261
     &        l_diag,k_diag
 
262
      integer g_tran, ! in
 
263
     &        g_anl,  ! in/out
 
264
     &        g_tmpanl! scratch GA array
 
265
      logical oprint
 
266
      character*(256) cstemp
 
267
      double precision sum
 
268
      double precision ga_trace_diag,trace,ac_trace
 
269
      external ga_trace_diag
 
270
        if (.not. ma_push_get(mt_dbl,nocc,'diag',l_diag,k_diag))
 
271
     &        call errquit('error alloc MA diag', 0, MA_ERR)
 
272
c Note.- This routine is for testing only
 
273
c ====>  it will be better to use the idea of 
 
274
c        trace() defined in aor_r1_tensor.F
 
275
c        test: symmetrize the g_anl matrix before the LMO trafo:
 
276
            if (oprint)
 
277
     &         write (luout,*) 'Message from beta_anl: Symmetrizing X'
 
278
            call ga_symmetrize(g_anl)
 
279
 
 
280
            write (cstemp,'(a)') 'aor_beta: tmpanl'
 
281
            if(.not.ga_create(MT_DBL,nocc,nocc,trim(cstemp),
 
282
     &         -1,-1,g_tmpanl))
 
283
     &         call errquit (trim(cstemp),0,GA_ERR)           
 
284
            call ga_zero(g_tmpanl)
 
285
            call ga_dgemm('t', 'n',nocc,nocc,nocc, 
 
286
     $                    1.0d0,g_tran,g_anl,0.0d0,
 
287
     &                    g_tmpanl) ! out: g_tmpanl
 
288
            call ga_zero(g_anl)
 
289
c ========= Added by FA to compute trace ====== START
 
290
         trace=0.0d0
 
291
         do i=1,nocc
 
292
          ac_trace=ga_ddot_patch(g_tmpanl,'n',i,i   ,1,nocc,
 
293
     &                           g_tran  ,'n',1,nocc,i,i) 
 
294
          trace=trace+ac_trace
 
295
 
 
296
             if (ga_nodeid().eq.0) then
 
297
               write(*,2) i,trace
 
298
 2             format('trace(',i3,')=',f15.8)
 
299
              endif
 
300
 
 
301
         enddo ! end-loop-i
 
302
c ========= Added by FA to compute trace ====== END
 
303
            call ga_dgemm('n', 'n',nocc,nocc,nocc, 
 
304
     $                    1.0d0,g_tmpanl,g_tran,0.0d0,
 
305
     &                    g_anl) ! out: g_anl
 
306
            if (.not.ga_destroy(g_tmpanl)) call errquit
 
307
     &         ('aor_beta: ga_destroy failed g_tmpanl',0,GA_ERR)
 
308
 
 
309
            if (oprint) write (luout,
 
310
     &         '(/t12,a,t26,a/t11,6(''-''),t22,12(''-''))')
 
311
     &         'LMO #','contrib.'
 
312
            
 
313
            call ga_get_diagonal(g_anl,dbl_mb(k_diag))
 
314
            
 
315
            sum = 0.0d0
 
316
            do i=1,nocc
 
317
 
 
318
              sum = sum + dbl_mb(k_diag+i-1)
 
319
 
 
320
              if (ga_nodeid().eq.0) then
 
321
               write(*,1) i,sum
 
322
 1             format('sum(',i3,')=',f15.8)
 
323
              endif
 
324
 
 
325
              if (oprint) write (luout,'(t11,i6,t22,f12.4)')
 
326
     &           i,dbl_mb(k_diag+i-1)
 
327
            enddo ! end-loop-i
 
328
 
 
329
          if (ga_nodeid().eq.0) then
 
330
           write(*,23) 1,sum,ga_trace_diag(g_anl),trace
 
331
 23        format('g_anl-trace-lmo(',i3,')=(',
 
332
     &            f15.8,',',f15.8,',',f15.8,')')
 
333
          endif
 
334
 
 
335
            if (oprint)
 
336
     &         write (luout,'(1x,a,2i1,a,f12.4)')
 
337
     &         'Component ',idir,idir,': Sum = ',sum
 
338
 
 
339
            if (oprint) write (luout,'(1x,40(''-''))')
 
340
            
 
341
              sum = ga_trace_diag(g_anl)
 
342
            if (oprint) write (luout,*) 'sum from ga_trace: ',sum
 
343
 
 
344
            if (.not. ma_pop_stack(l_diag))
 
345
     &       call errquit('error deloc MA diag',0, MA_ERR)                  
 
346
      return
 
347
      end
 
348
 
 
349
      subroutine get_tracelessQtensor(g_quadel, ! in/out: quadrupole tensor
 
350
     &                                nbf,      ! in: nr. basis functions
 
351
     &                                g_work)   ! in: scratch GA array
 
352
c
 
353
c Author : Fredy W. Aquino
 
354
c Date   : 03-15-12
 
355
c Note.- Modified from original aoresponse source code
 
356
c        for extension to spin-unrestricted case
 
357
c        original aoresponse source code was written by 
 
358
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
359
c --> Experimental (not published yet)
 
360
        
 
361
      implicit none
 
362
#include "errquit.fh"
 
363
#include "global.fh"
 
364
#include "mafdecls.fh"
 
365
#include "msgids.fh"
 
366
#include "stdio.fh"
 
367
      integer g_quadel,g_work
 
368
      integer nbf,iresp
 
369
      integer alo(3),ahi(3), 
 
370
     &        blo(3),bhi(3)
 
371
      integer qindex(3,3)
 
372
      double precision tenm8,one,two,three,
 
373
     &                 zero,half,third,four
 
374
      parameter (tenm8=1d-8, one=1d0, two=2d0, three=3d0,
 
375
     &           zero=0d0, half=one/two,
 
376
     &           third=one/three, four=4.0d0)
 
377
c     define translation table for quadrupole indices in
 
378
c     packed storage
 
379
c     XX=1, XY=YX=2, XZ=ZX=3, YY=4, YZ=ZY=5, ZZ=6
 
380
      qindex(1,1) = 1
 
381
      qindex(1,2) = 2
 
382
      qindex(2,1) = 2
 
383
      qindex(1,3) = 3
 
384
      qindex(3,1) = 3
 
385
      qindex(2,2) = 4
 
386
      qindex(2,3) = 5
 
387
      qindex(3,2) = 5
 
388
      qindex(3,3) = 6   
 
389
c     a) from trace -> g_work
 
390
      call ga_zero(g_work)  ! use for trace
 
391
      do iresp = 1,3          
 
392
        alo(1) = 1
 
393
        ahi(1) = nbf
 
394
        alo(2) = 1
 
395
        ahi(2) = nbf
 
396
        alo(3) = qindex(iresp,iresp) 
 
397
        ahi(3) = qindex(iresp,iresp)
 
398
        blo(1) = 1
 
399
        bhi(1) = nbf
 
400
        blo(2) = 1
 
401
        bhi(2) = nbf           
 
402
        call nga_add_patch(one,g_quadel,alo,ahi,
 
403
     &                     one,g_work  ,blo,bhi,
 
404
     &                         g_work  ,blo,bhi)
 
405
       enddo ! end-loop-iresp         
 
406
c      b) scale quadel by 3
 
407
       call ga_scale(g_quadel,three)         
 
408
c      c) subtract trace from diagonal        
 
409
       do iresp = 1,3           
 
410
        alo(1) = 1
 
411
        ahi(1) = nbf
 
412
        alo(2) = 1
 
413
        ahi(2) = nbf
 
414
        alo(3) = qindex(iresp,iresp) 
 
415
        ahi(3) = qindex(iresp,iresp)
 
416
        blo(1) = 1
 
417
        bhi(1) = nbf
 
418
        blo(2) = 1
 
419
        bhi(2) = nbf          
 
420
        call nga_add_patch(one,g_quadel,alo,ahi,
 
421
     &                    -one,g_work  ,blo,bhi,
 
422
     &                         g_quadel,alo,ahi)
 
423
       enddo ! end-loop-iresp         
 
424
c      d) divide the result by two, then by three
 
425
c         because of the factor by which the result enters
 
426
c         the Buckingham-Dunn tensor        
 
427
       call ga_scale(g_quadel,half*third)        
 
428
      return
 
429
      end
 
430
 
 
431
      subroutine write_vects(
 
432
     &              rtdb,geom,basis, ! in: handles
 
433
     &              g_vecE1,         ! in: C^(1,E)
 
434
     &              g_vecB1,         ! in: C^(1,B)
 
435
     &              g_tran,          ! in:
 
436
     &              g_vectmp,        ! in: scratch GA array
 
437
     &              npol,nocc,       ! in: nr. polariz,occ MOs
 
438
     &              nocct,nvirt,     ! in: nr. occ,virt MOs
 
439
     &              nopen,nmot,      ! in: nr. open shells,nmot=nocc*nvirt
 
440
     &              froct,epst,      ! in:
 
441
     &              nbf,             ! in: nr. basis functions
 
442
     &              lmo,             ! in: logical flag
 
443
     &              debug)           ! in: logical flag for debugging
 
444
c
 
445
c Author : Fredy W. Aquino
 
446
c Date   : 03-15-12
 
447
c Note.- Modified from original aoresponse source code
 
448
c        for extension to spin-unrestricted case
 
449
c        original aoresponse source code was written by 
 
450
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
451
c --> Experimental (not published yet)
 
452
        
 
453
      implicit none
 
454
#include "errquit.fh"
 
455
#include "global.fh"
 
456
#include "mafdecls.fh"
 
457
#include "msgids.fh"
 
458
#include "geom.fh"
 
459
#include "rtdb.fh"
 
460
#include "bas.fh"
 
461
#include "stdio.fh"
 
462
#include "prop.fh"
 
463
      integer rtdb,geom,basis
 
464
      integer g_vecE1,g_vecB1,g_tran,
 
465
     &        g_vectmp(2)
 
466
      integer idir,nbf,npol,nocc
 
467
      logical lmo,debug
 
468
      integer alo(3),ahi(3), 
 
469
     &        blo(3),bhi(3),
 
470
     &        clo(3),chi(3)
 
471
      character*(256) cstemp
 
472
      integer nocct(npol),nvirt(npol),
 
473
     &        nopen(npol),nmot(npol)
 
474
      double precision froct(nbf,2),epst(nbf,2)
 
475
c        ---------------------------------------------------
 
476
c        part of the analysis is storing the perturbed
 
477
c        MOs. As elsewhere, we assume a closed shell system.
 
478
c        We write a lot of data here, e-field and b-field
 
479
c        perturbed canonical MOs and, if applicable, LMOs
 
480
c       ----------------------------------------------------
 
481
        do idir = 1,3
 
482
          alo(1) = 1
 
483
          ahi(1) = nbf
 
484
          alo(2) = 1
 
485
          ahi(2) = nbf
 
486
          alo(3) = idir         
 
487
          ahi(3) = idir
 
488
          blo(1) = 1
 
489
          bhi(1) = nbf
 
490
          blo(2) = 1
 
491
          bhi(2) = nbf
 
492
          call ga_zero(g_vectmp(1))
 
493
          call nga_copy_patch('n',g_vecE1    ,alo,ahi,
 
494
     &                            g_vectmp(1),blo,bhi)
 
495
          if (debug) then
 
496
           if (ga_nodeid().eq.0) then
 
497
            write(*,24) idir
 
498
 24         format('---- g_vectmp-vecE1(',i3,')-------- START')
 
499
           endif
 
500
           call ga_print(g_vectmp(1))
 
501
           if (ga_nodeid().eq.0) then
 
502
            write(*,25) idir
 
503
 25         format('---- g_vectmp-vecE1(',i3,')-------- END')
 
504
           endif
 
505
          endif ! end-if-debug
 
506
          write(cstemp,'(a,i1,a)') 'cmo-efield',idir,'.movecs'
 
507
          call hnd_vec_write(rtdb,geom,basis,
 
508
     &                       nbf,nocct,nopen,nvirt,'scf',
 
509
     &                       g_vectmp,froct,epst,nmot,cstemp)
 
510
          call ga_zero(g_vectmp(1))
 
511
          call nga_copy_patch('n',g_vecB1    ,alo,ahi,
 
512
     &                            g_vectmp(1),blo,bhi)
 
513
          if (debug) then
 
514
           if (ga_nodeid().eq.0) then
 
515
            write(*,26) idir
 
516
 26         format('---- g_vectmp-vecB1(',i3,')-------- START')
 
517
           endif
 
518
           call ga_print(g_vectmp(1))
 
519
           if (ga_nodeid().eq.0) then
 
520
            write(*,27) idir
 
521
 27         format('---- g_vectmp-vecB1(',i3,')-------- END')
 
522
           endif
 
523
          endif ! end-if-debug
 
524
          write(cstemp,'(a,i1,a)') 'cmo-bfield',idir,'.movecs'
 
525
          call hnd_vec_write(rtdb,geom,basis,
 
526
     &                       nbf,nocct,nopen,nvirt,'scf',
 
527
     &                       g_vectmp,froct,epst,nmot,cstemp)
 
528
          if (lmo) then
 
529
            alo(1) = 1
 
530
            ahi(1) = nbf
 
531
            alo(2) = 1
 
532
            ahi(2) = nocc
 
533
            alo(3) = idir         
 
534
            ahi(3) = idir
 
535
            blo(1) = 1
 
536
            bhi(1) = nocc
 
537
            blo(2) = 1
 
538
            bhi(2) = nocc
 
539
            clo(1) = 1
 
540
            chi(1) = nbf
 
541
            clo(2) = 1
 
542
            chi(2) = nocc
 
543
            call ga_zero(g_vectmp(1))
 
544
            call nga_matmul_patch('n','n',1d0,0d0,
 
545
     &                            g_vecE1    ,alo,ahi,
 
546
     &                            g_tran     ,blo,bhi,
 
547
     &                            g_vectmp(1),clo,chi)
 
548
            if (debug) then
 
549
             if (ga_nodeid().eq.0) then
 
550
              write(*,28) idir
 
551
 28           format('---- g_vectmp-vecE1-lmo(',i3,')---- START')
 
552
             endif
 
553
             call ga_print(g_vectmp(1))
 
554
             if (ga_nodeid().eq.0) then
 
555
              write(*,29) idir
 
556
 29           format('---- g_vectmp-vecE1-lmo(',i3,')---- END')
 
557
             endif
 
558
            endif ! end-if-debug
 
559
            write(cstemp,'(a,i1,a)') 'lmo-efield',idir,'.movecs'
 
560
            call hnd_vec_write(rtdb,geom,basis,
 
561
     &                         nbf,nocct,nopen,nvirt,'scf',
 
562
     &                         g_vectmp,froct,epst,nmot,cstemp)
 
563
            call ga_zero(g_vectmp(1))
 
564
            call nga_matmul_patch('n','n',1d0,0d0,
 
565
     &                            g_vecB1    ,alo,ahi,
 
566
     &                            g_tran     ,blo,bhi,
 
567
     &                            g_vectmp(1),clo,chi)
 
568
            if (debug) then
 
569
             if (ga_nodeid().eq.0) then
 
570
              write(*,30) idir
 
571
 30           format('---- g_vectmp-vecB1-lmo(',i3,')---- START')
 
572
             endif
 
573
             call ga_print(g_vectmp(1))
 
574
             if (ga_nodeid().eq.0) then
 
575
              write(*,31) idir
 
576
 31           format('---- g_vectmp-vecB1-lmo(',i3,')---- END')
 
577
             endif
 
578
            endif ! end-if-debug
 
579
            write(cstemp,'(a,i1,a)') 'lmo-bfield',idir,'.movecs'
 
580
            call hnd_vec_write(rtdb,geom,basis,
 
581
     &                         nbf,nocct,nopen,nvirt,'scf',
 
582
     &                         g_vectmp,froct,epst,nmot, cstemp)
 
583
          end if ! lmo
 
584
        enddo ! idir
 
585
      return
 
586
      end