← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/property/aor_r1_beta_anl.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/property/aor_r1_beta_anl.F
 
1
 
 
2
      subroutine aor_r1_beta_anl (rtdb, basis, geom, omega, lstatic,
 
3
     &   ncomp,g_smat0, g_sket1, g_vecB1,  g_dipel, g_quadel, g_vectors,
 
4
     &   froct, epst, nbf, nmo, nocct, nvirt, lgiao, lquad, lanalyze,
 
5
     &   lvelocity ,lifetime, lmagpert, g_vecE1, g_vecE1_im)
 
6
 
 
7
c $Id: aor_r1_beta_anl.F 21176 2011-10-10 06:35:49Z d3y133 $
 
8
      
 
9
c     =================================================================
 
10
      
 
11
c     purpose: analyze beta tensor. See routine
 
12
c     aor_r1_beta.F for additional comments. 
 
13
c     Use a molecular orientation in which the chiral
 
14
c     response tensor is diagonal!
 
15
 
 
16
c     called from: aoresponse_driver_new
 
17
 
 
18
c     =================================================================
 
19
 
 
20
      implicit none
 
21
 
 
22
#include "errquit.fh"
 
23
#include "global.fh"
 
24
#include "mafdecls.fh"
 
25
#include "msgids.fh"
 
26
#include "geom.fh"
 
27
#include "rtdb.fh"
 
28
#include "bas.fh"
 
29
#include "stdio.fh"
 
30
#include "apiP.fh"
 
31
#include "prop.fh"
 
32
#include "bgj.fh"
 
33
 
 
34
c     ---------------------
 
35
c     subroutine arguments:
 
36
c     ---------------------
 
37
 
 
38
      integer rtdb    ! [input] run-time database handle
 
39
      integer basis   ! [input] basis handle
 
40
      integer geom    ! [input] geometry handle
 
41
 
 
42
c     These are all input, too
 
43
      integer g_smat0, g_vectors(2), g_dipel, 
 
44
     &   g_sket1, g_quadel, g_vecB1, g_vecE1(2), g_vecE1_im(2)
 
45
      integer nfreq, response_order, nbf, nmo, ncomp
 
46
      integer nocct(2), nvirt(2)
 
47
      double precision froct(nbf,2), epst(nbf,2)
 
48
      double precision gamwidth, omega
 
49
      logical lgiao, lquad, lanalyze, lvelocity, lifetime, lmagpert,
 
50
     &   lstatic
 
51
 
 
52
c     ----------------
 
53
c     local variables:
 
54
c     ----------------
 
55
 
 
56
c     global array handles:
 
57
      
 
58
      integer
 
59
     &   g_anl, g_work,
 
60
     &   g_temp, g_tmpanl, g_tran, g_vectmp(2)
 
61
 
 
62
      integer l_diag, k_diag
 
63
 
 
64
c     other local variables: 
 
65
 
 
66
      integer nmot(2), nocvir(2), nopen(2)
 
67
      data nopen(1),nopen(2)/0,0/
 
68
 
 
69
      integer dims(3), chunk(3)
 
70
      integer alo(3), ahi(3), blo(3), bhi(3), clo(3), chi(3)
 
71
 
 
72
c     dipole-quadrupole polarizability, cartesian rep.:
 
73
      double precision dipquadre(3,6), dipquadim(3,6)
 
74
c     traceless dipole-quadrupole tensor, full storage
 
75
      double precision dqpol(3,3,3)
 
76
 
 
77
      integer LCTensor(3,3,3)
 
78
      integer qindex(3,3)
 
79
      double precision tmpmat(3,3)
 
80
 
 
81
      character*(256) cstemp
 
82
 
 
83
      character*(1) direction(3)
 
84
      data direction/'x','y','z'/
 
85
      
 
86
      integer ispin, nspin
 
87
      integer ipm, nocc, nvir, nocv, imo, jmo, nmo1, idir, iresp,
 
88
     &   i,j,k,l
 
89
      logical debug, dbgmat, 
 
90
     &   lzora, lantisym, lmo, status, oprint
 
91
      double precision sum, scaling
 
92
      double precision tenm8, one, two, three, zero, half, third, four
 
93
      parameter (tenm8=1d-8, one=1d0, two=2d0, three=3d0,
 
94
     &   zero=0d0, half=one/two,
 
95
     &   third=one/three, four=4.0d0)
 
96
 
 
97
c     external functions:
 
98
 
 
99
      character*(256) lmotrans
 
100
      logical file_read_ga
 
101
      external file_read_ga
 
102
 
 
103
      double precision ga_trace_diag
 
104
      external ga_trace_diag
 
105
 
 
106
c  ====================================================================
 
107
 
 
108
      debug = .false.  ! .true. during development
 
109
      dbgmat = .false. ! debug large matrices
 
110
      oprint = ga_nodeid().eq.0
 
111
 
 
112
      if (debug) write (luout,*) 'hello from aor_r1_beta_anl'
 
113
 
 
114
      dipquadre(:,:) = 0
 
115
      dipquadim(:,:) = 0
 
116
 
 
117
c     make sure lvelocity.ne.T., we do not support that in this
 
118
c     subroutine to keep the clutter at a manageable level.
 
119
c     same for lmagpert
 
120
 
 
121
      if (lvelocity) 
 
122
     &   call errquit ('aor_beta: lvelocity set',1,INPUT_ERR)
 
123
 
 
124
      if (lmagpert) 
 
125
     &   call errquit ('aor_beta: lmagpert set',1,INPUT_ERR)
 
126
 
 
127
 
 
128
c     -------------------------
 
129
c     define Levi-Civita tensor for quadrupole additions
 
130
c     -------------------------
 
131
      LCtensor(:,:,:) = 0      
 
132
      LCtensor(1,2,3) = 1
 
133
      LCtensor(2,3,1) = 1
 
134
      LCtensor(3,1,2) = 1      
 
135
      LCtensor(2,1,3) = -1
 
136
      LCtensor(1,3,2) = -1
 
137
      LCtensor(3,2,1) = -1        
 
138
 
 
139
c     define translation table for quadrupole incices in
 
140
c     packed storage
 
141
c     XX=1, XY=YX=2, XZ=ZX=3, YY=4, YZ=ZY=5, ZZ=6
 
142
      qindex(1,1) = 1
 
143
      qindex(1,2) = 2
 
144
      qindex(2,1) = 2
 
145
      qindex(1,3) = 3
 
146
      qindex(3,1) = 3
 
147
      qindex(2,2) = 4
 
148
      qindex(2,3) = 5
 
149
      qindex(3,2) = 5
 
150
      qindex(3,3) = 6      
 
151
 
 
152
c     set parameters that control the various computational options
 
153
c     (later we will set most of this by input)
 
154
      nspin      =  1           ! assume closed shell
 
155
      lzora      = .false.      ! not yet available here 
 
156
 
 
157
      if (debug) write (luout,*) 'giao, velocity',
 
158
     &    lgiao, lvelocity
 
159
 
 
160
 
 
161
c     -----------------------------------------
 
162
c     determine number of occ * virt orbitals
 
163
c     and nmot(1:2) and fix froct, if necessary
 
164
c     -----------------------------------------
 
165
 
 
166
      do ispin = 1,nspin
 
167
        nocvir(ispin) = nocct(ispin) * nvirt(ispin)
 
168
        nmot(ispin) = nmo
 
169
        if (nmo .lt.nbf) then
 
170
          do imo = nmo+1,nbf
 
171
            froct(imo,ispin) = 0d0
 
172
          enddo
 
173
        endif
 
174
      end do
 
175
 
 
176
c     ----------------------------------------------
 
177
c     check if we have localized MOs on file. If yes
 
178
c     read them, assuming nspin.eq.1
 
179
c     ----------------------------------------------
 
180
 
 
181
      i = 1
 
182
      status=rtdb_get(rtdb,'prop:pmlocalization',MT_INT,1,i)
 
183
      lmo = (i.eq.0)
 
184
      if (lmo) then
 
185
        if (oprint) write (luout,*) 'analysis: LMO switch found'
 
186
        write (cstemp,'(a)') 'aor_beta: g_tran'
 
187
        if(.not.ga_create(MT_DBL, nocct(1), nocct(1), trim(cstemp),
 
188
     &     -1,-1,g_tran))
 
189
     &     call errquit (trim(cstemp),0,GA_ERR)
 
190
        if (debug) write (luout,*) 'g_tran allocated'
 
191
        call util_file_name('lmotrans',.true.,.true.,lmotrans)
 
192
        if(.not.file_read_ga(lmotrans,g_tran)) call errquit
 
193
     $     ('aor_r1_beta_anl: could not read lmotrans',0, DISK_ERR)
 
194
      end if ! lmo
 
195
 
 
196
c     ----------------------------------------------------
 
197
c     start loop over spins (nspin=2 not yet supported !!!)
 
198
c     ----------------------------------------------------
 
199
      
 
200
      do ispin = 1, nspin
 
201
        
 
202
        nmo1 = nmot(ispin)      ! total no.of MOs for this spin
 
203
        nocc = nocct(ispin)     ! occupied MOs
 
204
        nvir = nvirt(ispin)     ! virtual MOs
 
205
        nocv = nocvir(ispin)    ! nocc * nvir
 
206
        
 
207
c       ------------------------------
 
208
c       allocate some temp. work space
 
209
c       ------------------------------
 
210
        
 
211
        chunk(1) = nbf
 
212
        chunk(2) = -1
 
213
        dims(1) = nbf
 
214
        dims(2) = nbf        
 
215
        write(cstemp,'(a)') 'work'
 
216
        if (.not.nga_create(MT_DBL,2,dims,cstemp(1:4),chunk,
 
217
     &     g_work)) call 
 
218
     &     errquit('aoresponse: nga_create failed: '//cstemp(1:4),
 
219
     &     0,GA_ERR)     
 
220
        call ga_zero (g_work)
 
221
        
 
222
c       allocate intermediate vector for matrix multiplications
 
223
c       used to create the final results
 
224
        
 
225
        write (cstemp,'(a)') 'aor_beta: temp1'
 
226
        if(.not.ga_create(MT_DBL, nbf, nocc, trim(cstemp),
 
227
     &     -1,-1,g_temp))
 
228
     &     call errquit (trim(cstemp),0,GA_ERR)
 
229
        if (debug) write (luout,*) 'g_temp allocated'
 
230
 
 
231
c       allocate matrix that accumulates the analysis data
 
232
 
 
233
        write (cstemp,'(a)') 'aor_beta: g_anl'
 
234
        if(.not.ga_create(MT_DBL, nocc, nocc, trim(cstemp),
 
235
     &     -1,-1,g_anl))
 
236
     &     call errquit (trim(cstemp),0,GA_ERR)
 
237
        if (debug) write (luout,*) 'g_anl allocated'
 
238
        call ga_zero(g_anl)
 
239
 
 
240
c       diagonal elements of the last matrix
 
241
        if (.not. ma_push_get(mt_dbl, nocc, 'diag', l_diag, k_diag))
 
242
     &        call errquit('error alloc MA diag', 0, MA_ERR)
 
243
 
 
244
c       lmos: debug transformation
 
245
        if (lmo .and. debug) then
 
246
          call ga_print(g_tran)
 
247
          call ga_dgemm('t', 'n', nocc, nocc, nocc, 
 
248
     $      1.0d0, g_tran, g_tran, 0.0d0, g_anl)
 
249
          call ga_print(g_anl)
 
250
        end if
 
251
        
 
252
c       ---------------------------------------------------------
 
253
c       solution of CPKS is in g_vecE1. Below we need
 
254
c       only the sum of the +/- components so we add them here
 
255
c       and store them in g_vecE1(1)
 
256
c       ---------------------------------------------------------
 
257
        
 
258
        if (ncomp.gt.1) then
 
259
          call ga_add(1d0, g_vecE1(1), 1d0,  g_vecE1(2),
 
260
     &       g_vecE1(1))
 
261
          if (lifetime) then
 
262
            call ga_add(1d0, g_vecE1_im(1), 1d0,  g_vecE1_im(2),
 
263
     &         g_vecE1_im(1))
 
264
          end if
 
265
        endif    
 
266
 
 
267
c       ------------------------------------------------
 
268
c       for Buckingham-Dunn tensor we need the traceless
 
269
c       quadrupole tensor
 
270
c       ------------------------------------------------
 
271
 
 
272
        if (lquad) then
 
273
          
 
274
c         a) from trace -> g_work
 
275
          call ga_zero(g_work)  ! use for trace
 
276
          do iresp = 1,3
 
277
            
 
278
            alo(1) = 1
 
279
            ahi(1) = nbf
 
280
            alo(2) = 1
 
281
            ahi(2) = nbf
 
282
            alo(3) = qindex(iresp,iresp) 
 
283
            ahi(3) = qindex(iresp,iresp)
 
284
            blo(1) = 1
 
285
            bhi(1) = nbf
 
286
            blo(2) = 1
 
287
            bhi(2) = nbf
 
288
            
 
289
            call nga_add_patch(one,g_quadel,alo,ahi,
 
290
     &         one,g_work,blo,bhi,g_work,blo,bhi)
 
291
          end do                ! iresp
 
292
          
 
293
c         b) scale quadel by 3
 
294
          call ga_scale(g_quadel,three)
 
295
          
 
296
c         c) subtract trace from diagonal
 
297
          
 
298
          do iresp = 1,3
 
299
            
 
300
            alo(1) = 1
 
301
            ahi(1) = nbf
 
302
            alo(2) = 1
 
303
            ahi(2) = nbf
 
304
            alo(3) = qindex(iresp,iresp) 
 
305
            ahi(3) = qindex(iresp,iresp)
 
306
            blo(1) = 1
 
307
            bhi(1) = nbf
 
308
            blo(2) = 1
 
309
            bhi(2) = nbf
 
310
            
 
311
            call nga_add_patch(one,g_quadel,alo,ahi,
 
312
     &         -one,g_work,blo,bhi,g_quadel,alo,ahi)
 
313
          end do                ! iresp
 
314
          
 
315
c         d) divide the result by two, then by three
 
316
c         because of the factor by which the result enters
 
317
c         the Buckingham-Dunn tensor
 
318
          
 
319
          call ga_scale(g_quadel, half*third)
 
320
          
 
321
!          if (debug) call ga_print(g_quadel)
 
322
        end if                  ! lquad
 
323
        
 
324
c       ---------------------------------------------------------
 
325
c       start loop over the components of the response tensor and 
 
326
c       calculate the final results
 
327
c       ---------------------------------------------------------
 
328
        
 
329
        do idir = 1,3           ! direction of the perturbing field
 
330
 
 
331
c         g_anl is going to accumulate the results
 
332
          call ga_zero(g_anl)
 
333
 
 
334
          if (oprint)
 
335
     &       write (luout,'(1x,40(''-'')/1x,a,2i1)')
 
336
     &       'MO analysis of OR tensor component ',idir,idir
 
337
                                   
 
338
c           -------------------------------------------------------
 
339
c           (A) calculate optical rotation beta from C(E) S(0) C(B)
 
340
c           ------------------------------------------------------
 
341
            
 
342
            alo(1) = 1
 
343
            ahi(1) = nbf
 
344
            alo(2) = 1
 
345
            ahi(2) = nbf
 
346
            alo(3) = 1 
 
347
            ahi(3) = 1
 
348
            blo(1) = 1
 
349
            bhi(1) = nbf
 
350
            blo(2) = 1
 
351
            bhi(2) = nocc
 
352
            blo(3) = idir      ! pick magnetic field direction
 
353
            bhi(3) = idir 
 
354
            clo(1) = 1
 
355
            chi(1) = nbf
 
356
            clo(2) = 1
 
357
            chi(2) = nocc
 
358
            
 
359
            call ga_zero(g_temp)
 
360
            call nga_matmul_patch('n','n',1d0,0d0,
 
361
     &         g_smat0,alo,ahi,
 
362
     &         g_vecB1,blo,bhi,
 
363
     &         g_temp,clo,chi)
 
364
            
 
365
            if (debug) write (luout,*)
 
366
     &         'beta: S(0) C(B) intermediate complete'
 
367
            
 
368
            alo(1) = 1
 
369
            ahi(1) = nocc
 
370
            alo(2) = 1
 
371
            ahi(2) = nbf
 
372
            alo(3) = idir
 
373
            ahi(3) = idir
 
374
            blo(1) = 1
 
375
            bhi(1) = nbf
 
376
            blo(2) = 1
 
377
            bhi(2) = nocc
 
378
            clo(1) = 1
 
379
            chi(1) = nocc
 
380
            clo(2) = 1
 
381
            chi(2) = nocc
 
382
            
 
383
            call ga_zero(g_work)
 
384
            call nga_matmul_patch('t','n',1d0,0d0,
 
385
     &         g_vecE1,alo,ahi,
 
386
     &         g_temp,blo,bhi,
 
387
     &         g_work,clo,chi)
 
388
                        
 
389
c           the factor of two is for the orbital occupations,
 
390
c           assuming that ispin is never equal to two
 
391
            
 
392
            scaling = two
 
393
            if (lstatic) scaling = four
 
394
            call nga_add_patch(scaling,g_work,clo,chi,
 
395
     &         one,g_anl,clo,chi,g_anl,clo,chi)
 
396
 
 
397
            if (debug) write (luout,*) 'beta: C(E) S(0) C(B) complete'
 
398
            
 
399
            if (lgiao) then
 
400
              
 
401
c             --------------------------------------
 
402
c             if we use GIAOs there is a second term
 
403
c             in beta which is C(E) S(1ket) C(0)
 
404
c             --------------------------------------
 
405
              
 
406
              alo(1) = 1
 
407
              ahi(1) = nbf
 
408
              alo(2) = 1
 
409
              ahi(2) = nbf
 
410
              alo(3) = idir    ! pick the correct sket1 direction
 
411
              ahi(3) = idir
 
412
              blo(1) = 1
 
413
              bhi(1) = nbf
 
414
              blo(2) = 1
 
415
              bhi(2) = nocc
 
416
              clo(1) = 1
 
417
              chi(1) = nbf
 
418
              clo(2) = 1
 
419
              chi(2) = nocc
 
420
                        
 
421
              call ga_zero(g_temp)
 
422
              call nga_matmul_patch('n','n',1d0,0d0,
 
423
     &           g_sket1,alo,ahi,
 
424
     &           g_vectors(ispin),blo,bhi,
 
425
     &           g_temp,clo,chi)
 
426
              
 
427
              if (debug) write (luout,*)
 
428
     &           'beta: S(ket1) C(0) intermediate complete'
 
429
                            
 
430
              alo(1) = 1
 
431
              ahi(1) = nocc
 
432
              alo(2) = 1
 
433
              ahi(2) = nbf
 
434
              alo(3) = idir
 
435
              ahi(3) = idir
 
436
              blo(1) = 1
 
437
              bhi(1) = nbf
 
438
              blo(2) = 1
 
439
              bhi(2) = nocc
 
440
              clo(1) = 1
 
441
              chi(1) = nocc
 
442
              clo(2) = 1
 
443
              chi(2) = nocc
 
444
              
 
445
              call ga_zero(g_work)
 
446
              call nga_matmul_patch('t','n',1d0,0d0,
 
447
     &           g_vecE1,alo,ahi,
 
448
     &           g_temp,blo,bhi,
 
449
     &           g_work,clo,chi)
 
450
                            
 
451
c             the factor of two is for the orbital occupations,
 
452
c             assuming that ispin is never equal to two
 
453
 
 
454
              scaling = two
 
455
              if (lstatic) scaling = four
 
456
              call nga_add_patch(scaling,g_work,clo,chi,
 
457
     &           one,g_anl,clo,chi,g_anl,clo,chi)
 
458
              
 
459
              if (debug) write (luout,*)
 
460
     &           'beta: C(E) S(ket1) C(0) complete'
 
461
 
 
462
            end if              ! lgiao
 
463
            
 
464
 
 
465
 
 
466
c           ----------------------------------------------------
 
467
c           if requested by input, add to OR beta the quadrupole
 
468
c           polarizability terms
 
469
c           ----------------------------------------------------       
 
470
            
 
471
            if (lquad) then
 
472
              
 
473
c             add the quadrupole terms to the B.-D. tensor
 
474
              do k = 1,3
 
475
                do l = 1,3
 
476
                  
 
477
                  if (LCtensor(idir,k,l).eq.0) goto 1000
 
478
                  
 
479
                  iresp = qindex(l,idir)
 
480
                  
 
481
                  alo(1) = 1
 
482
                  ahi(1) = nbf
 
483
                  alo(2) = 1
 
484
                  ahi(2) = nbf
 
485
                  alo(3) = iresp ! pick direction iresp for g_quadel
 
486
                  ahi(3) = iresp
 
487
                  blo(1) = 1
 
488
                  bhi(1) = nbf
 
489
                  blo(2) = 1
 
490
                  bhi(2) = nocc
 
491
                  clo(1) = 1
 
492
                  chi(1) = nbf
 
493
                  clo(2) = 1
 
494
                  chi(2) = nocc
 
495
                
 
496
                
 
497
                  call ga_zero(g_temp)
 
498
                  call nga_matmul_patch('n','n',1d0,0d0,
 
499
     &               g_quadel,alo,ahi,
 
500
     &               g_vectors(ispin),blo,bhi,
 
501
     &               g_temp,clo,chi)
 
502
                  
 
503
                  if (debug) write (luout,*)
 
504
     &               'quad: h(Q) C(0) intermediate complete'
 
505
                                
 
506
                  alo(1) = 1
 
507
                  ahi(1) = nocc
 
508
                  alo(2) = 1
 
509
                  ahi(2) = nbf
 
510
                  alo(3) = k ! not idir, see equation
 
511
                  ahi(3) = k
 
512
                  blo(1) = 1
 
513
                  bhi(1) = nbf
 
514
                  blo(2) = 1
 
515
                  bhi(2) = nocc
 
516
                  clo(1) = 1
 
517
                  chi(1) = nocc
 
518
                  clo(2) = 1
 
519
                  chi(2) = nocc
 
520
                  
 
521
                  call ga_zero(g_work)
 
522
                  
 
523
                  call nga_matmul_patch('t','n',1d0,0d0,
 
524
     &               g_vecE1,alo,ahi,
 
525
     &               g_temp,blo,bhi,
 
526
     &               g_work,clo,chi)
 
527
                  
 
528
c                 the factor of two is for the orbital occupations,
 
529
c                 assuming that ispin is never equal to two
 
530
                  
 
531
                  scaling = -two
 
532
                  if (lstatic) scaling = -four
 
533
                  
 
534
c                 Levi-Civita symbol:
 
535
                  scaling = scaling * LCtensor(idir,k,l)
 
536
                  
 
537
                  if (debug) write (luout,*) 'scaling=',scaling
 
538
                  call nga_add_patch(scaling,g_work,clo,chi,
 
539
     &               one,g_anl,clo,chi,g_anl,clo,chi)
 
540
                
 
541
                if (debug) write (luout,*)
 
542
     &             'quad C(Q) h(E) C(0) complete'
 
543
 
 
544
 1000           continue ! jump here if LCTensor(idir,k,l) is 0
 
545
                
 
546
              end do            ! l
 
547
            end do              ! k
 
548
 
 
549
          end if                ! lquad
 
550
 
 
551
c         -----------------------------------------
 
552
c         end loop over responding field components
 
553
c         -----------------------------------------
 
554
 
 
555
c         ---------------------
 
556
c         Canonical MO analysis
 
557
c         ---------------------
 
558
 
 
559
          if (oprint) write (luout,
 
560
     &       '(/t12,a,t26,a/t11,6(''-''),t22,12(''-''))')
 
561
     &       'CMO #','contrib.'
 
562
 
 
563
          call ga_get_diagonal(g_anl, dbl_mb(k_diag) )
 
564
 
 
565
          sum = zero
 
566
          do i = 1,nocc
 
567
            sum = sum + dbl_mb(k_diag+i-1)
 
568
            if (oprint) write (luout,'(t11,i6,t22,f12.4)')
 
569
     &         i,dbl_mb(k_diag+i-1)
 
570
          end do                ! i
 
571
          if (oprint)
 
572
     &       write (luout,'(1x,a,2i1,a,f12.4)') 'Component ',idir,idir,
 
573
     &       ': Sum = ',sum
 
574
 
 
575
          if (oprint) write (luout,'(1x,40(''-''))')
 
576
 
 
577
          if (debug) then
 
578
            sum = ga_trace_diag(g_anl)
 
579
            if (oprint) write (luout,*) 'sum from ga_trace: ',sum
 
580
          end if                ! debug
 
581
 
 
582
c         ---------------------
 
583
c         Localized MO analysis
 
584
c         ---------------------
 
585
 
 
586
          if (lmo) then
 
587
 
 
588
c           test: symmetrize the g_anl matrix before the LMO trafo:
 
589
            if (oprint)
 
590
     &         write (luout,*) 'Message from beta_anl: Symmetrizing X'
 
591
            call ga_symmetrize(g_anl)
 
592
 
 
593
            write (cstemp,'(a)') 'aor_beta: tmpanl'
 
594
            if(.not.ga_create(MT_DBL, nocc, nocc, trim(cstemp),
 
595
     &         -1,-1,g_tmpanl))
 
596
     &         call errquit (trim(cstemp),0,GA_ERR)
 
597
            if (debug) write (luout,*) 'g_tmpanl allocated'
 
598
            
 
599
            call ga_zero(g_tmpanl)
 
600
            call ga_dgemm('t', 'n', nocc, nocc, nocc, 
 
601
     $         1.0d0, g_tran, g_anl, 0.0d0, g_tmpanl)
 
602
            call ga_zero(g_anl)
 
603
            call ga_dgemm('n', 'n', nocc, nocc, nocc, 
 
604
     $         1.0d0, g_tmpanl, g_tran, 0.0d0, g_anl)
 
605
 
 
606
          if (.not.ga_destroy(g_tmpanl))
 
607
     &         call errquit
 
608
     &         ('aor_beta: ga_destroy failed g_tmpanl',
 
609
     &         0,GA_ERR)
 
610
 
 
611
            if (oprint) write (luout,
 
612
     &         '(/t12,a,t26,a/t11,6(''-''),t22,12(''-''))')
 
613
     &         'LMO #','contrib.'
 
614
            
 
615
            call ga_get_diagonal(g_anl, dbl_mb(k_diag) )
 
616
            
 
617
            sum = zero
 
618
            do i = 1,nocc
 
619
              sum = sum + dbl_mb(k_diag+i-1)
 
620
              if (oprint) write (luout,'(t11,i6,t22,f12.4)')
 
621
     &           i,dbl_mb(k_diag+i-1)
 
622
            end do              ! i
 
623
            if (oprint)
 
624
     &         write (luout,'(1x,a,2i1,a,f12.4)')
 
625
     &         'Component ',idir,idir,': Sum = ',sum
 
626
 
 
627
            if (oprint) write (luout,'(1x,40(''-''))')
 
628
            
 
629
            if (debug) then
 
630
              sum = ga_trace_diag(g_anl)
 
631
              if (oprint) write (luout,*) 'sum from ga_trace: ',sum
 
632
            end if              ! debug
 
633
            
 
634
            
 
635
          end if                ! lmo
 
636
          
 
637
        end do                  ! idir = 1,3
 
638
 
 
639
c       -------------------------------------------
 
640
c       end loop over perturbing E-field components
 
641
c       -------------------------------------------
 
642
 
 
643
 
 
644
 
 
645
c       -----------------
 
646
c       deallocate memory
 
647
c       -----------------
 
648
 
 
649
          if (.not.ga_destroy(g_temp))
 
650
     &       call errquit
 
651
     &       ('aor_beta: ga_destroy failed g_temp',
 
652
     &       0,GA_ERR)
 
653
 
 
654
 
 
655
        if (.not.ga_destroy(g_work))
 
656
     &     call 
 
657
     &     errquit('aoresponse: ga_destroy failed g_work',
 
658
     &     0,GA_ERR)
 
659
 
 
660
          if (.not.ga_destroy(g_anl))
 
661
     &       call errquit
 
662
     &       ('aor_beta: ga_destroy failed g_anl',
 
663
     &       0,GA_ERR)
 
664
 
 
665
         if (.not. ma_pop_stack(l_diag))
 
666
     &       call errquit('error deloc MA diag',0, MA_ERR)
 
667
 
 
668
c        ---------------------------------------------------
 
669
c        part of the analysis is storing the perturbed
 
670
c        MOs. As elsewhere, we assume a closed shell system.
 
671
c        We write a lot of data here, e-field and b-field
 
672
c        perturbed canonical MOs and, if applicable, LMOs
 
673
c       ----------------------------------------------------
 
674
 
 
675
         chunk(1) = nbf
 
676
         chunk(2) = -1
 
677
         dims(1) = nbf
 
678
         dims(2) = nbf        
 
679
         write(cstemp,'(a)') 'vectmp(1)'
 
680
         if (.not.nga_create(MT_DBL,2,dims,cstemp(1:4),chunk,
 
681
     &      g_vectmp(1))) call 
 
682
     &     errquit('aoresponse: nga_create failed: '//cstemp(1:9),
 
683
     &     0,GA_ERR)     
 
684
        call ga_zero (g_vectmp(1))
 
685
 
 
686
        do idir = 1,3
 
687
 
 
688
          alo(1) = 1
 
689
          ahi(1) = nbf
 
690
          alo(2) = 1
 
691
          ahi(2) = nbf
 
692
          alo(3) = idir         
 
693
          ahi(3) = idir
 
694
          blo(1) = 1
 
695
          bhi(1) = nbf
 
696
          blo(2) = 1
 
697
          bhi(2) = nbf
 
698
 
 
699
 
 
700
          call ga_zero(g_vectmp(1))
 
701
          call nga_copy_patch('n',g_vecE1,alo,ahi,g_vectmp(1),blo,bhi)
 
702
 
 
703
          write(cstemp,'(a,i1,a)') 'cmo-efield',idir,'.movecs'
 
704
          call hnd_vec_write(rtdb,geom,basis,nbf,nocct,nopen, nvirt
 
705
     &       ,'scf',g_vectmp,froct, epst,nmot, cstemp)
 
706
 
 
707
          call ga_zero(g_vectmp(1))
 
708
          call nga_copy_patch('n',g_vecB1,alo,ahi,g_vectmp(1),blo,bhi)
 
709
 
 
710
          write(cstemp,'(a,i1,a)') 'cmo-bfield',idir,'.movecs'
 
711
          call hnd_vec_write(rtdb,geom,basis,nbf,nocct,nopen, nvirt
 
712
     &       ,'scf',g_vectmp,froct, epst,nmot, cstemp)
 
713
 
 
714
          if (lmo) then
 
715
 
 
716
            alo(1) = 1
 
717
            ahi(1) = nbf
 
718
            alo(2) = 1
 
719
            ahi(2) = nocct(1)
 
720
            alo(3) = idir         
 
721
            ahi(3) = idir
 
722
            blo(1) = 1
 
723
            bhi(1) = nocct(1)
 
724
            blo(2) = 1
 
725
            bhi(2) = nocct(1)
 
726
            clo(1) = 1
 
727
            chi(1) = nbf
 
728
            clo(2) = 1
 
729
            chi(2) = nocct(1)
 
730
 
 
731
            call ga_zero(g_vectmp(1))
 
732
            call nga_matmul_patch('n','n',1d0,0d0,
 
733
     &         g_vecE1,alo,ahi,
 
734
     &         g_tran,blo,bhi,
 
735
     &         g_vectmp(1),clo,chi)
 
736
 
 
737
            write(cstemp,'(a,i1,a)') 'lmo-efield',idir,'.movecs'
 
738
            call hnd_vec_write(rtdb,geom,basis,nbf,nocct,nopen, nvirt
 
739
     &         ,'scf',g_vectmp,froct, epst,nmot, cstemp)
 
740
 
 
741
            call ga_zero(g_vectmp(1))
 
742
            call nga_matmul_patch('n','n',1d0,0d0,
 
743
     &         g_vecB1,alo,ahi,
 
744
     &         g_tran,blo,bhi,
 
745
     &         g_vectmp(1),clo,chi)
 
746
 
 
747
            write(cstemp,'(a,i1,a)') 'lmo-bfield',idir,'.movecs'
 
748
            call hnd_vec_write(rtdb,geom,basis,nbf,nocct,nopen, nvirt
 
749
     &         ,'scf',g_vectmp,froct, epst,nmot, cstemp)
 
750
 
 
751
          end if ! lmo
 
752
 
 
753
        end do ! idir
 
754
 
 
755
        if (.not.ga_destroy(g_vectmp(1))) call errquit ('aor_beta:
 
756
     &     ga_destroy failed g_vectmp(1)', 0,GA_ERR)
 
757
 
 
758
c       -------------------------------------------------
 
759
c       un-add the frequency components in vec_E1 in case
 
760
c       we reuse these arrays:
 
761
c       -------------------------------------------------
 
762
 
 
763
        if (ncomp.gt.1) then
 
764
          call ga_add(1d0, g_vecE1(1), -1d0,  g_vecE1(2),
 
765
     &       g_vecE1(1))
 
766
          if (lifetime) then
 
767
            call ga_add(1d0, g_vecE1_im(1), -1d0,  g_vecE1_im(2),
 
768
     &         g_vecE1_im(1))
 
769
          end if
 
770
        endif     
 
771
 
 
772
        
 
773
      enddo                     ! ispin = 1,2 from way above
 
774
              
 
775
c     ---------------------------------------------------------------
 
776
c     end loop over spin components (which we don't use right now
 
777
c     since nspin is forced to be 1 at the beginning of this routine)
 
778
c     ---------------------------------------------------------------
 
779
            
 
780
 
 
781
      if (lmo) then
 
782
        if (.not.ga_destroy(g_tran))
 
783
     &       call errquit
 
784
     &       ('aor_beta: ga_destroy failed g_tran',0,GA_ERR)
 
785
      end if
 
786
 
 
787
c     ----------------
 
788
c     all done. return
 
789
c     ----------------
 
790
                  
 
791
      
 
792
c     ==================================================================
 
793
      
 
794
      return
 
795
      
 
796
      end
 
797