← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/property/CalcPerturbedTDPmat1_opt.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/CalcPerturbedTDPmat1_opt.F
 
1
      subroutine CalcPerturbedTDPmat1_opt(
 
2
     &                 ncomp,    ! in : nr. components to calculate
 
3
     &                 g_pmats,  ! out: density matrix      symmetrized
 
4
     &                 g_pmata,  ! out: density matrix  antisymmetrized -NOT USED
 
5
     &                 g_amat,   ! in :   perturbed MO coefficients 
 
6
     &                 g_vectors,! in : unperturbed MO coefficients
 
7
     &                 naos,     ! in : nr. AOs
 
8
     &                 nocc,     ! in : nr. occupied MOs
 
9
     &                 nvir,     ! in : nr. virtual  MOs
 
10
     &                 nmo,      ! in : nr. MOs
 
11
     &                 lantisym, ! in : = .true. calc. (symm,antisymm)=(pmats,pmata)
 
12
     &                 lstatic,  ! in : = .true. static response, dynamic otherwise
 
13
     &                 imag,     ! in : = .true. if amat is imaginary instead of real
 
14
     &                 haveocc)  ! in : = .true. if amat contains occ-occ block
 
15
 
 
16
* $Id: CalcPerturbedTDPmat1.F 19707 2010-10-29 17:59:36Z d3y133 $
 
17
 
 
18
c     ==================================================================
 
19
 
 
20
c     calculate frequency-dependent density matrix perturbation
 
21
c     (symmetric and antisymmetric part), linear response,
 
22
c     from a set of perturbed MO coefficients
 
23
 
 
24
c     we assume DOUBLE occupation of all occupied orbitals and
 
25
c     REAL unperturbed orbitals. The perturbation can be either
 
26
c     purely real or purely imaginary
 
27
 
 
28
c     THIS ROUTINE USES TOO MUCH MEMORY; IT COULD DO HE SAME
 
29
C     JOB WITH LESS TEMP SPACE. FIX THIS
 
30
 
 
31
c     input: 
 
32
c     ncomp      - number of components to be calculated
 
33
c     g_amat     - the perturbed MO coefficients are
 
34
c                  written as C(+/-) = C(0)*A(+/-),
 
35
c                  g_amat contains the elements of matrix A
 
36
c                  (only the virt - occ block, or nmo - occ block)
 
37
c     g_vectors  - unperturbed MO coefficients C(0)
 
38
c     lantisym   - logical switch to calculate symmetric
 
39
c                  and antisymmetric
 
40
c                  part separately or just the total density matrix
 
41
c     lstatic    - static response, assume that both components
 
42
c                  of amat are equal. assumes ncomp = 1 (!)
 
43
c     imag      - true if amat is imaginary instad of real
 
44
c     haveocc    - true if amat contains occ-occ block, too
 
45
      
 
46
c     output : g_pmats, g_pmata: symmetric and antisymmetric
 
47
c     part of perturbed density matrix, global arrays, if (lantisym),
 
48
c     otherwise the total density matrix is in pmats, and pmata=0
 
49
      
 
50
c     remark : all perturbed matrices are classified by
 
51
c     (+/-) frequency components 
 
52
c     
 
53
c     remark: the density matrix is given by
 
54
c     transpose(P) = C n C(dagger), i.e. in the end we transpose the
 
55
c     result to get the correct density matrix out
 
56
 
 
57
c     ==================================================================
 
58
c
 
59
c Author : Fredy W. Aquino
 
60
c Date   : 03-15-12
 
61
c Note.- Modified from original aoresponse source code
 
62
c        for extension to spin-unrestricted case
 
63
c        original aoresponse source code was written by 
 
64
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
65
c        Modifying original source code: CalcPerturbedTDPmat1()
 
66
c --> Experimental (not published yet)
 
67
    
 
68
      implicit none
 
69
#include "errquit.fh"
 
70
#include "global.fh"
 
71
#include "mafdecls.fh"
 
72
c     subroutine arguments:
 
73
      integer ncomp
 
74
      integer g_pmats(ncomp), 
 
75
     &        g_pmata(ncomp), 
 
76
     &        g_amat(ncomp),
 
77
     &        g_vectors          
 
78
      integer naos, nocc, nvir, nmo
 
79
      logical lantisym
 
80
      logical lstatic, imag, haveocc
 
81
 
 
82
c     local variables:
 
83
      integer ip, im, ipm
 
84
      double precision half, one, two,toscl
 
85
      parameter (half = 0.5d0, one = 1.0d0)
 
86
      logical debug
 
87
      external get_PertDens
 
88
 
 
89
c      debug = .true. !     allow debugging printout
 
90
      debug = .false. ! not allow debugging printout
 
91
c     check range of ncomp
 
92
      if (ncomp.le.0 .or. ncomp.gt.2) then
 
93
        call errquit('CalcPerturbedTDPmat: ncomp out of range',
 
94
     &     0,CALC_ERR)
 
95
      endif
 
96
      if (ncomp.gt.1 .and. lstatic) then
 
97
        call errquit
 
98
     &     ('CalcPerturbedTDPmat1: ncomp > 1 but lstatic.eq.true.',
 
99
     &     0,CALC_ERR)
 
100
      endif
 
101
c     assign + and - components for indexing amat:
 
102
      if (lstatic) then
 
103
        ip = 1
 
104
        im = 1
 
105
      else
 
106
        ip = 1
 
107
        im = 2
 
108
      endif
 
109
 
 
110
      call ga_sync()
 
111
 
 
112
c      if (ga_nodeid().eq.0) then
 
113
c       write(*,190) ip,im,lstatic,imag
 
114
c 190     format('(ip,im,lstatic,imag)=(',i3,',',i3,',',L1,',',L1,')')
 
115
c      endif
 
116
        
 
117
      do ipm = 1,ncomp
 
118
        call ga_zero(g_pmats(ipm))
 
119
c        call ga_zero(g_pmata(ipm))
 
120
      enddo
 
121
 
 
122
      if (nocc+nvir .ne. nmo) call errquit
 
123
     &   ('CalcPerturbedTDPmat1: wrong no. of orbitals',0,CALC_ERR)
 
124
 
 
125
c     -------------------------------------------------------------
 
126
c     First we assemble P(+). Note that A(-) is assumed to be A(-)*
 
127
c     in fact (A = amat)
 
128
c     This allows us to use the same algorithm no matter if A is
 
129
c     real and symmetric or imaginary and antisymmetric
 
130
c     -------------------------------------------------------------
 
131
c     ----------------------------
 
132
c     First  step: C    n C(-,dagger)
 
133
c     Second step: C(+) n C(dagger)
 
134
c     -----------------------------
 
135
      if (debug) then
 
136
       do ipm=1,2
 
137
        if (ga_nodeid().eq.0)
 
138
     &   write(*,*) '---- BEF get_PertDens g_amat(',ipm,')--- START'
 
139
         call ga_print(g_amat(ipm))
 
140
        if (ga_nodeid().eq.0)
 
141
     &   write(*,*) '---- BEF get_PertDens g_amat(',ipm,')--- END'
 
142
       enddo ! end-loop-ipm 
 
143
      endif ! end-if-debug     
 
144
 
 
145
      call get_PertDens(
 
146
     &               g_pmats(1),!out:      symmetrized dens
 
147
     &               g_pmata(1),!out: anti symmetrized dens
 
148
     &               g_amat,    ! in: the u mat
 
149
     &               g_vectors, ! in: MO vect
 
150
     &               ncomp,     ! in: nr. components g_amat
 
151
     &               ip,im,     ! in: indices of u vect
 
152
     &               imag,      ! in: = T -> imag
 
153
     &               haveocc,   ! in: logical flag
 
154
     &               lantisym,  ! in: logical flag
 
155
     &               naos,nmo,  ! in: nr. AOs,MOs
 
156
     &               nocc,nvir, ! in: nr. (occ,virt) MOs
 
157
     &               debug)     ! in: =.true. -> show debug printouts
 
158
 
 
159
        if (debug) then
 
160
             if (ga_nodeid().eq.0)
 
161
     &       write(*,*) '---- g_pmats-nw-------- START'
 
162
             call ga_print(g_pmats(1))
 
163
            if (ga_nodeid().eq.0)
 
164
     &       write(*,*) '---- g_pmats-nw-------- END'      
 
165
c             if (ga_nodeid().eq.0)
 
166
c     &       write(*,*) '---- g_pmata-nw-------- START'
 
167
c             call ga_print(g_pmata(1))
 
168
c            if (ga_nodeid().eq.0)
 
169
c     &       write(*,*) '---- g_pmata-nw-------- END'   
 
170
        endif ! end-if-debug     
 
171
      if (lstatic .or. ncomp.eq.1) then
 
172
c       skip calculation of component 2 of the density matrix
 
173
        if (ga_nodeid().eq.0)
 
174
     &   write(*,*) 'FA-Skipping calc of 2nd component'
 
175
 
 
176
        goto 7000
 
177
      endif  
 
178
c        if (ga_nodeid().eq.0)
 
179
c     &   write(*,*) 'FA-calculating 2nd component'
 
180
c     ----------------------------
 
181
c     First  step: C    n C(+,dagger)
 
182
c     Second step: C(-) n C(dagger)
 
183
c     -----------------------------
 
184
c     Note 1.- I tried calling get_PertDens 2nd time
 
185
c            and swapping (ip,im) --> (im,ip)
 
186
c            and I got the relationship shown below:
 
187
c     lantisym=.true. :
 
188
c     g_pmats(2)= -transpose(g_pmats(1))
 
189
c     g_pmata(2)= g_pmata(1)= 0
 
190
c     lantisym=.false. :
 
191
c     g_pmats(2)= -transpose(g_pmats(1))
 
192
c     g_pmata(2)= g_pmata(1)
 
193
c     Note 2.- Swapping (ip,im) --> (im,ip)
 
194
c              is equivalent to doing the second
 
195
c              part (e.g. calc. of P(-)
 
196
      call ga_transpose(g_pmats(1),g_pmats(2))
 
197
                toscl= 1.0d0
 
198
      if (imag) toscl=-1.0d0
 
199
c      call ga_scale(g_pmats(2),-1.0d0)
 
200
      call ga_scale(g_pmats(2),toscl)
 
201
c      call ga_copy(g_pmata(1),g_pmata(2))
 
202
c     jump here from above in case of static calculation
 
203
 7000 continue    
 
204
      return
 
205
      end
 
206
 
 
207
      subroutine get_PertDens(
 
208
     &               g_pmats,  !out:      symmetrized dens
 
209
     &               g_pmata,  !out: anti symmetrized dens -NOT USED
 
210
     &               g_amat,   ! in: the u mat
 
211
     &               g_vectors,! in: MO vect
 
212
     &               ncomp,    ! in: nr. components g_amat
 
213
     &               ip,im,    ! in: indices of u vect
 
214
     &               imag,     ! in: = T -> imag
 
215
     &               haveocc,  ! in: logical flag
 
216
     &               lantisym, ! in: logical flag
 
217
     &               naos,nmo, ! in: nr. AOs,MOs
 
218
     &               nocc,nvir,! in: nr. (occ,virt) MOs
 
219
     &               debug)    ! in: =.true. -> show debug printouts
 
220
c
 
221
c Author : Fredy W. Aquino
 
222
c Date   : 03-15-12
 
223
c Note.- Modified from original aoresponse source code
 
224
c        for extension to spin-unrestricted case
 
225
c        original aoresponse source code was written by 
 
226
c        J. Autschbach and appears on nwchem-devtrunk (date:03-02-12)
 
227
c --> Experimental (not published yet)
 
228
 
 
229
      implicit none
 
230
#include "errquit.fh"
 
231
#include "global.fh"
 
232
#include "mafdecls.fh"
 
233
 
 
234
c     subroutine arguments:
 
235
c Note.- For (im,ip) --> (g_pmats,g_pmata)(1)
 
236
c        For (ip,im) --> (g_pmats,g_pmata)(2)
 
237
      integer ncomp,ip,im
 
238
      integer g_pmats,g_pmata,
 
239
     &        g_amat(ncomp),
 
240
     &        g_vectors           
 
241
      integer naos,nocc,nvir,nmo
 
242
      logical lantisym,imag,haveocc
 
243
c     local variables:
 
244
      integer g_ptmp,g_eig1,g_work
 
245
      double precision half, one, two
 
246
      parameter (half = 0.5d0, one = 1.0d0)
 
247
      logical debug
 
248
c     ------------------------
 
249
c     allocate workspace (GAs)
 
250
c     ------------------------    
 
251
      if (.not. ga_create(
 
252
     &    MT_DBL,naos,naos,'get_PertDens:ptmp',0,0,g_ptmp)) 
 
253
     &    call errquit('get_PertDens:ptmp',0,GA_ERR)
 
254
      if (.not. ga_create(
 
255
     &    MT_DBL,naos,naos,'get_PertDens:work',0,0,g_work)) 
 
256
     &    call errquit('get_PertDens:work',0,GA_ERR)
 
257
      if (.not. ga_create(
 
258
     &    MT_DBL,naos,naos,'get_PertDens:eig1',0,0,g_eig1)) 
 
259
     &    call errquit('get_PertDens:eig1',0,GA_ERR)
 
260
c     ----------------------------
 
261
c     First step:  C n C(-,dagger)
 
262
c     ----------------------------
 
263
c     calculate C(-,dagger)
 
264
      two = 2d0
 
265
      if (imag) two = -2d0
 
266
      if (.not.haveocc) then
 
267
 
 
268
       if (debug) then
 
269
         if (ga_nodeid().eq.0)
 
270
     &   write(*,*) 'FA-enter-no-haveocc...'
 
271
         if (ga_nodeid().eq.0)
 
272
     &    write(*,*) '-----g_vectors-nohaveocc---- START'
 
273
          call ga_print(g_vectors)
 
274
         if (ga_nodeid().eq.0)
 
275
     &    write(*,*) '-----g_vectors-nohaveocc---- END'
 
276
         if (ga_nodeid().eq.0) then
 
277
           write(*,1) im
 
278
 1         format('-----g_amat(',i3,')-nohaveocc---- START')
 
279
         endif
 
280
         call ga_print(g_amat(im))
 
281
         if (ga_nodeid().eq.0) then
 
282
           write(*,2) im
 
283
 2         format('-----g_amat(',i3,')-nohaveocc---- END')
 
284
         endif
 
285
          if (ga_nodeid().eq.0) then
 
286
           write(*,35) ip
 
287
 35         format('-----g_amat(',i3,')-yeshaveocc---- START')
 
288
          endif
 
289
          call ga_print(g_amat(ip))
 
290
          if (ga_nodeid().eq.0) then
 
291
           write(*,36) ip
 
292
 36         format('-----g_amat(',i3,')-yeshaveocc---- END')
 
293
          endif
 
294
 
 
295
       endif ! end-if-debug
 
296
 
 
297
       call ga_zero(g_eig1)
 
298
       call ga_matmul_patch('n','n', two,0d0,
 
299
     &                      g_vectors ,1,naos,nocc+1,nmo,
 
300
     &                      g_amat(im),1,nvir,1     ,nocc,
 
301
     &                      g_eig1    ,1,naos,1     ,nocc)
 
302
        if (debug) then
 
303
         if (ga_nodeid().eq.0)
 
304
     &    write(*,*) '-----g_eig1-nohaveocc---- START'
 
305
          call ga_print(g_eig1)
 
306
         if (ga_nodeid().eq.0)
 
307
     &    write(*,*) '-----g_eig1-nohaveocc---- END'
 
308
        endif ! end-if-debug
 
309
      else
 
310
        if (debug) then
 
311
          if (ga_nodeid().eq.0)
 
312
     &    write(*,*) '-----g_vectors-yeshaveocc---- START'
 
313
          call ga_print(g_vectors)
 
314
         if (ga_nodeid().eq.0)
 
315
     &    write(*,*) '-----g_vectors-yeshaveocc---- END'
 
316
          if (ga_nodeid().eq.0) then
 
317
           write(*,3) im
 
318
 3         format('-----g_amat(',i3,')-yeshaveocc---- START')
 
319
          endif
 
320
          call ga_print(g_amat(im))
 
321
          if (ga_nodeid().eq.0) then
 
322
           write(*,4) im
 
323
 4         format('-----g_amat(',i3,')-yeshaveocc---- END')
 
324
          endif
 
325
          if (ga_nodeid().eq.0) then
 
326
           write(*,33) ip
 
327
 33         format('-----g_amat(',i3,')-yeshaveocc---- START')
 
328
          endif
 
329
          call ga_print(g_amat(ip))
 
330
          if (ga_nodeid().eq.0) then
 
331
           write(*,34) ip
 
332
 34         format('-----g_amat(',i3,')-yeshaveocc---- END')
 
333
          endif
 
334
 
 
335
        endif ! end-if-debug
 
336
 
 
337
       call ga_zero(g_eig1)
 
338
       call ga_matmul_patch('n','n', two,0d0,
 
339
     &                      g_vectors ,1,naos,1,nmo,
 
340
     &                      g_amat(im),1,nmo ,1,nocc,
 
341
     &                      g_eig1    ,1,naos,1,nocc)
 
342
 
 
343
         if (debug) then
 
344
          if (ga_nodeid().eq.0)
 
345
     &     write(*,*) '-----g_eig1-yeshaveocc---- START'
 
346
          call ga_print(g_eig1)
 
347
          if (ga_nodeid().eq.0)
 
348
     &     write(*,*) '-----g_eig1-yeshaveocc---- END'
 
349
         endif ! end-if-debug
 
350
      endif ! .not.haveocc
 
351
c     note: the dimensioning for array B is that of the transposed
 
352
c     matrix, not of the original matrix. 
 
353
      
 
354
c     calculate C(0)C(-,dagger), store in g_ptmp
 
355
 
 
356
      call ga_zero(g_ptmp)
 
357
      call ga_matmul_patch('n','t', 1d0,0d0,
 
358
     &                     g_vectors,1,naos,1,nocc,
 
359
     &                     g_eig1   ,1,nocc,1,naos,
 
360
     &                     g_ptmp   ,1,naos,1,naos)   
 
361
  
 
362
c     -----------------------------
 
363
c     Second step: C(+) n C(dagger)
 
364
c     -----------------------------     
 
365
c     calculate C(+)
 
366
      two = 2.0d0
 
367
      if (.not.haveocc) then
 
368
 
 
369
       call ga_zero(g_eig1)
 
370
       call ga_matmul_patch('n','n', two,0d0,
 
371
     &                      g_vectors ,1,naos,nocc+1,nmo,
 
372
     &                      g_amat(ip),1,nvir,1     ,nocc,
 
373
     &                      g_eig1    ,1,naos,1     ,nocc)
 
374
 
 
375
       if (debug) then
 
376
        if (ga_nodeid().eq.0)
 
377
     &   write(*,*) '-----g_eig1-nohaveocc-2-- START'
 
378
        call ga_print(g_eig1)
 
379
        if (ga_nodeid().eq.0)
 
380
     &   write(*,*) '-----g_eig1-nohaveocc-2-- END'
 
381
       endif
 
382
      else
 
383
 
 
384
       call ga_zero(g_eig1)
 
385
       call ga_matmul_patch('n','n', two,0d0,
 
386
     &                      g_vectors ,1,naos,1,nmo,
 
387
     &                      g_amat(ip),1,nmo ,1,nocc,
 
388
     &                      g_eig1    ,1,naos,1,nocc)
 
389
 
 
390
       if (debug) then
 
391
        if (ga_nodeid().eq.0)
 
392
     &   write(*,*) '-----g_eig1-yeshaveocc-2---- START'
 
393
        call ga_print(g_eig1)
 
394
        if (ga_nodeid().eq.0)
 
395
     &   write(*,*) '-----g_eig1-yeshaveocc-2---- END'
 
396
       endif ! end-if-debug
 
397
      endif ! end-if-haveocc
 
398
      
 
399
c     calculate C(+)C(0,dagger), store in g_work
 
400
 
 
401
      call ga_zero(g_work)
 
402
      call ga_matmul_patch('n','t', 1d0,0d0,
 
403
     &                     g_eig1   ,1,naos,1,nocc,
 
404
     &                     g_vectors,1,nocc,1,naos,
 
405
     &                     g_work   ,1,naos,1,naos)
 
406
 
 
407
c     add the two terms together and transpose the density matrix
 
408
      if (debug) then
 
409
             if (ga_nodeid().eq.0)
 
410
     &       write(*,*) '---- CC1^t-------- START'
 
411
             call ga_print(g_ptmp)
 
412
            if (ga_nodeid().eq.0)
 
413
     &       write(*,*) '---- CC1^t--------- END'           
 
414
             if (ga_nodeid().eq.0)
 
415
     &       write(*,*) '---- C1C^t-------- START'
 
416
             call ga_print(g_work)
 
417
            if (ga_nodeid().eq.0)
 
418
     &       write(*,*) '---- C1C^t--------- END'     
 
419
      endif ! end-if-debug
 
420
c Doing: C(Du_q)^T + (Du_p)C^T
 
421
c     (q,p)=(im,ip) --> P(+)
 
422
c     (q,p)=(ip,im) --> P(-)
 
423
      call ga_add(1d0,g_ptmp,1d0,g_work,g_work) 
 
424
      call ga_transpose(g_work, g_ptmp)
 
425
c     calculate symmetrized and antisymmetrized part (+ component)
 
426
c     if requested on input:
 
427
      if (lantisym) then
 
428
        call ga_transpose(g_ptmp,g_work)
 
429
        call ga_add(half,g_ptmp, half,g_work,g_pmats)
 
430
      else
 
431
        call ga_copy(g_ptmp, g_pmats)
 
432
      endif
 
433
      if (debug) then
 
434
             if (ga_nodeid().eq.0)
 
435
     &       write(*,*) '---- g_pmats-------- START'
 
436
             call ga_print(g_pmats)
 
437
            if (ga_nodeid().eq.0)
 
438
     &       write(*,*) '---- g_pmats-------- END'    
 
439
c             if (ga_nodeid().eq.0)
 
440
c     &       write(*,*) '---- g_pmata-------- START'
 
441
c             call ga_print(g_pmata)
 
442
c            if (ga_nodeid().eq.0)
 
443
c     &       write(*,*) '---- g_pmata-------- END'           
 
444
      endif ! end-if-debug  
 
445
      if (.not.ga_destroy(g_ptmp))
 
446
     &   call errquit('get_PertDens: ga_destroy failed g_ptmp',
 
447
     &                0,GA_ERR)
 
448
      if (.not.ga_destroy(g_work))
 
449
     &   call errquit('get_PertDens: ga_destroy failed g_work',
 
450
     &                0,GA_ERR)
 
451
      if (.not.ga_destroy(g_eig1))
 
452
     &   call errquit('get_PertDens: ga_destroy failed g_eig1',
 
453
     &                0,GA_ERR) 
 
454
 
 
455
      return
 
456
      end    
 
457
 
 
458
c ++++++++++++++++++ Added for debugging ++++++++++++++++ START
 
459
 
 
460
      subroutine CalcPerturbedTDPmat1_deb
 
461
     &   (ncomp, g_pmats, g_pmata, g_amat, g_vectors, naos, nocc,
 
462
     &   nvir, nmo, lantisym, lstatic, imag, haveocc)
 
463
 
 
464
* $Id: CalcPerturbedTDPmat1.F 19707 2010-10-29 17:59:36Z d3y133 $
 
465
 
 
466
c     ==================================================================
 
467
 
 
468
c     calculate frequency-dependent density matrix perturbation
 
469
c     (symmetric and antisymmetric part), linear response,
 
470
c     from a set of perturbed MO coefficients
 
471
 
 
472
c     we assume DOUBLE occupation of all occupied orbitals and
 
473
c     REAL unperturbed orbitals. The perturbation can be either
 
474
c     purely real or purely imaginary
 
475
 
 
476
c     THIS ROUTINE USES TOO MUCH MEMORY; IT COULD DO HE SAME
 
477
C     JOB WITH LESS TEMP SPACE. FIX THIS
 
478
 
 
479
c     input: 
 
480
c     ncomp      - number of components to be calculated
 
481
c     g_amat     - the perturbed MO coefficients are
 
482
c                  written as C(+/-) = C(0)*A(+/-),
 
483
c                  g_amat contains the elements of matrix A
 
484
c                  (only the virt - occ block, or nmo - occ block)
 
485
c     g_vectors  - unperturbed MO coefficients C(0)
 
486
c     lantisym   - logical switch to calculate symmetric
 
487
c                  and antisymmetric
 
488
c                  part separately or just the total density matrix
 
489
c     lstatic    - static response, assume that both components
 
490
c                  of amat are equal. assumes ncomp = 1 (!)
 
491
c     imag      - true if amat is imaginary instad of real
 
492
c     haveocc    - true if amat contains occ-occ block, too
 
493
      
 
494
c     output : g_pmats, g_pmata: symmetric and antisymmetric
 
495
c     part of perturbed density matrix, global arrays, if (lantisym),
 
496
c     otherwise the total density matrix is in pmats, and pmata=0
 
497
      
 
498
c     remark : all perturbed matrices are classified by
 
499
c     (+/-) frequency components 
 
500
c     
 
501
c     remark: the density matrix is given by
 
502
c     transpose(P) = C n C(dagger), i.e. in the end we transpose the
 
503
c     result to get the correct density matrix out
 
504
 
 
505
c     ==================================================================
 
506
      
 
507
      implicit none
 
508
 
 
509
#include "errquit.fh"
 
510
#include "global.fh"
 
511
#include "mafdecls.fh"
 
512
 
 
513
c     subroutine arguments:
 
514
      integer ncomp
 
515
      integer g_pmats(ncomp), g_pmata(ncomp), g_amat(*),
 
516
     &   g_vectors           ! GA
 
517
      integer naos, nocc, nvir, nmo
 
518
      logical lantisym
 
519
      logical lstatic, imag, haveocc
 
520
 
 
521
c     local variables:
 
522
      integer g_ptmp, g_eig1, g_work
 
523
      integer imo, imu, inu, ll, ip, im, ipm
 
524
      double precision rtemp
 
525
      double precision half, one, two
 
526
      parameter (half = 0.5d0, one = 1.0d0)
 
527
      logical debug
 
528
 
 
529
      integer type, dim1, dim2
 
530
      external get_PertDens
 
531
 
 
532
c     ==================================================================
 
533
c ++++++++++++ WARNING -CHECKING lantisym=T +++++++ START
 
534
c       if (ga_nodeid().eq.0)
 
535
c     &  write(*,*) 'FA-WARNING: check lantisym=T'
 
536
c       lantisym=.true.
 
537
c ++++++++++++ WARNING -CHECKING lantisym=T +++++++ TRUE
 
538
 
 
539
      debug = .false.
 
540
 
 
541
c     check range of ncomp
 
542
 
 
543
      if (ncomp.le.0 .or. ncomp.gt.2) then
 
544
        call errquit('CalcPerturbedTDPmat: ncomp out of range',
 
545
     &     0,CALC_ERR)
 
546
      endif
 
547
 
 
548
c     cowardy refuse so calculate two components of perturbed
 
549
c     density matrix if lstatic switch is set to true
 
550
 
 
551
      if (ncomp.gt.1 .and. lstatic) then
 
552
        call errquit
 
553
     &     ('CalcPerturbedTDPmat1: ncomp > 1 but lstatic.eq.true.',
 
554
     &     0,CALC_ERR)
 
555
      endif
 
556
 
 
557
 
 
558
c     assign + and - components for indexing amat:
 
559
 
 
560
      if (lstatic) then
 
561
        ip = 1
 
562
        im = 1
 
563
      else
 
564
        ip = 1
 
565
        im = 2
 
566
      endif
 
567
 
 
568
      if (ga_nodeid().eq.0) then
 
569
       write(*,190) ip,im,lstatic,imag
 
570
 190     format('(ip,im,lstatic,imag)=(',i3,',',i3,',',L1,',',L1,')')
 
571
      endif
 
572
        
 
573
      do ipm = 1,ncomp
 
574
        call ga_zero(g_pmats(ipm))
 
575
        call ga_zero(g_pmata(ipm))
 
576
      enddo
 
577
 
 
578
      if (debug) write (6,'(a,4i6)') 'nocc,nvir,nmo',nocc, nvir, nmo
 
579
 
 
580
      if (nocc+nvir .ne. nmo) call errquit
 
581
     &   ('CalcPerturbedTDPmat1: wrong no. of orbitals',0,CALC_ERR)
 
582
 
 
583
c     ------------------------
 
584
c     allocate workspace (GAs)
 
585
c     ------------------------
 
586
      
 
587
      if (.not. ga_create(MT_DBL, naos, naos,
 
588
     &   'CalcPerturbedTDPmat1:ptmp',
 
589
     &   0, 0, g_ptmp)) call errquit('CalcPerturbedTDPmat1:ptmp', 0,
 
590
     &   GA_ERR)
 
591
 
 
592
      if (.not. ga_create(MT_DBL, naos, naos,
 
593
     &   'CalcPerturbedTDPmat1:work',
 
594
     &   0, 0, g_work)) call errquit('CalcPerturbedTDPmat1:work', 0,
 
595
     &   GA_ERR)
 
596
      
 
597
      if (.not. ga_create(MT_DBL, naos, nocc,
 
598
     &   'CalcPerturbedTDPmat1:eig1',
 
599
     &   0, 0, g_eig1)) call errquit('CalcPerturbedTDPmat1:eig1', 0,
 
600
     &   GA_ERR)
 
601
 
 
602
      if (debug) then
 
603
c       debug array dimensions
 
604
        call ga_inquire (g_eig1,type, dim1, dim2)
 
605
        write (6,'(a,2i4)') 'g_eig1:',dim1,dim2
 
606
        call ga_inquire (g_ptmp,type, dim1, dim2)
 
607
        write (6,'(a,2i4)') 'g_ptmp:',dim1,dim2
 
608
        call ga_inquire (g_work,type, dim1, dim2)
 
609
        write (6,'(a,2i4)') 'g_work:',dim1,dim2
 
610
        call ga_inquire (g_amat(1),type, dim1, dim2)
 
611
        write (6,'(a,2i4)') 'g_amat(1):',dim1,dim2
 
612
        call ga_inquire (g_vectors,type, dim1, dim2)
 
613
        write (6,'(a,2i4)') 'g_vectors:',dim1,dim2
 
614
      endif
 
615
 
 
616
 
 
617
c     -------------------------------------------------------------
 
618
c     First we assemble P(+). Note that A(-) is assumed to be A(-)*
 
619
c     in fact (A = amat)
 
620
c     This allows us to use the same algorithm no matter if A is
 
621
c     real and symmetric or imaginary and antisymmetric
 
622
c     -------------------------------------------------------------
 
623
 
 
624
      call ga_zero(g_ptmp)
 
625
      call ga_zero(g_work)
 
626
      call ga_zero(g_eig1)
 
627
      call ga_sync()
 
628
 
 
629
c     ----------------------------
 
630
c     First step:  C n C(-,dagger)
 
631
c     ----------------------------
 
632
 
 
633
c     calculate C(-,dagger)
 
634
      two = 2d0
 
635
      if (imag) two = -2d0
 
636
 
 
637
       if (ga_nodeid().eq.0)
 
638
     &   write(*,31) imag,two
 
639
 31       format('FA-debug (imag,two)=(',L1,',',f15.8)
 
640
 
 
641
      if (.not.haveocc) then
 
642
        if (ga_nodeid().eq.0)
 
643
     &   write(*,*) 'FA-enter-no-haveocc...'
 
644
         if (ga_nodeid().eq.0)
 
645
     &    write(*,*) '-----g_vectors-nohaveocc---- START'
 
646
          call ga_print(g_vectors)
 
647
         if (ga_nodeid().eq.0)
 
648
     &    write(*,*) '-----g_vectors-nohaveocc---- END'
 
649
          if (ga_nodeid().eq.0) then
 
650
           write(*,1) im
 
651
 1         format('-----g_amat(',i3,')-nohaveocc---- START')
 
652
          endif
 
653
          call ga_print(g_amat(im))
 
654
          if (ga_nodeid().eq.0) then
 
655
           write(*,2) im
 
656
 2         format('-----g_amat(',i3,')-nohaveocc---- END')
 
657
          endif
 
658
          if (ga_nodeid().eq.0) then
 
659
           write(*,1070) ip
 
660
 1070         format('-----g_amat(',i3,')-nohaveocc---- START')
 
661
          endif
 
662
          call ga_print(g_amat(ip))
 
663
          if (ga_nodeid().eq.0) then
 
664
           write(*,1071) ip
 
665
 1071        format('-----g_amat(',i3,')-nohaveocc---- END')
 
666
          endif
 
667
       call ga_matmul_patch('n','n', two,0d0,
 
668
     &   g_vectors ,1,naos,nocc+1,nmo,
 
669
     &   g_amat(im),1,nvir,1     ,nocc,
 
670
     &   g_eig1    ,1,naos,1     ,nocc)
 
671
 
 
672
         if (ga_nodeid().eq.0)
 
673
     &    write(*,*) '-----g_eig1-nohaveocc---- START'
 
674
          call ga_print(g_eig1)
 
675
         if (ga_nodeid().eq.0)
 
676
     &    write(*,*) '-----g_eig1-nohaveocc---- END'
 
677
 
 
678
      else
 
679
         if (ga_nodeid().eq.0)
 
680
     &    write(*,*) 'FA-enter-yes-haveocc...'
 
681
 
 
682
         if (ga_nodeid().eq.0)
 
683
     &    write(*,*) '-----g_vectors-yeshaveocc---- START'
 
684
          call ga_print(g_vectors)
 
685
         if (ga_nodeid().eq.0)
 
686
     &    write(*,*) '-----g_vectors-yeshaveocc---- END'
 
687
          if (ga_nodeid().eq.0) then
 
688
           write(*,3) im
 
689
 3         format('-----g_amat(',i3,')-yeshaveocc---- START')
 
690
          endif
 
691
          call ga_print(g_amat(im))
 
692
          if (ga_nodeid().eq.0) then
 
693
           write(*,4) im
 
694
 4         format('-----g_amat(',i3,')-yeshaveocc---- END')
 
695
          endif
 
696
 
 
697
       call ga_matmul_patch('n','n', two,0d0,
 
698
     &   g_vectors ,1,naos,1,nmo,
 
699
     &   g_amat(im),1,nmo ,1,nocc,
 
700
     &   g_eig1    ,1,naos,1,nocc)
 
701
 
 
702
         if (ga_nodeid().eq.0)
 
703
     &    write(*,*) '-----g_eig1-yeshaveocc---- START'
 
704
         call ga_print(g_eig1)
 
705
         if (ga_nodeid().eq.0)
 
706
     &    write(*,*) '-----g_eig1-yeshaveocc---- END'
 
707
      endif ! .not.haveocc
 
708
      call ga_sync()
 
709
 
 
710
      if (debug) write (6,*) '1'
 
711
 
 
712
c     note: the dimensioning for array B is that of the transposed
 
713
c     matrix, not of the original matrix. 
 
714
      
 
715
c     calculate C(0)C(-,dagger), store in g_ptmp
 
716
      call ga_matmul_patch('n','t', 1d0,0d0,
 
717
     &   g_vectors,1,naos,1,nocc,
 
718
     &   g_eig1,1,nocc,1,naos,
 
719
     &   g_ptmp,1,naos,1,naos)  
 
720
      call ga_sync()
 
721
 
 
722
      if (debug) write (6,*) '2'   
 
723
      
 
724
c     -----------------------------
 
725
c     Second step: C(+) n C(dagger)
 
726
c     -----------------------------
 
727
      
 
728
c     calculate C(+)
 
729
      two = 2.0d0
 
730
      if (.not.haveocc) then
 
731
 
 
732
        if (ga_nodeid().eq.0)
 
733
     &   write(*,*) 'FA-enter-no-haveocc-2...'
 
734
 
 
735
       call ga_matmul_patch('n','n', two,0d0,
 
736
     &   g_vectors ,1,naos,nocc+1,nmo,
 
737
     &   g_amat(ip),1,nvir,1     ,nocc,
 
738
     &   g_eig1    ,1,naos,1     ,nocc)
 
739
 
 
740
       if (ga_nodeid().eq.0)
 
741
     &  write(*,*) '-----g_eig1-nohaveocc-2-- START'
 
742
       call ga_print(g_eig1)
 
743
       if (ga_nodeid().eq.0)
 
744
     &  write(*,*) '-----g_eig1-nohaveocc-2-- END'
 
745
 
 
746
      else
 
747
 
 
748
         if (ga_nodeid().eq.0)
 
749
     &   write(*,*) 'FA-enter-yes-haveocc-2...'
 
750
 
 
751
       call ga_matmul_patch('n','n', two,0d0,
 
752
     &   g_vectors ,1,naos,1,nmo,
 
753
     &   g_amat(ip),1,nmo ,1,nocc,
 
754
     &   g_eig1    ,1,naos,1,nocc)
 
755
 
 
756
        if (ga_nodeid().eq.0)
 
757
     &   write(*,*) '-----g_eig1-yeshaveocc-2---- START'
 
758
        call ga_print(g_eig1)
 
759
        if (ga_nodeid().eq.0)
 
760
     &   write(*,*) '-----g_eig1-yeshaveocc-2---- END'
 
761
 
 
762
      endif
 
763
      call ga_sync()
 
764
 
 
765
      if (debug) write (6,*) '3'
 
766
      
 
767
c     calculate C(+)C(0,dagger), store in g_work
 
768
      call ga_matmul_patch('n','t', 1d0,0d0,
 
769
     &   g_eig1   ,1,naos,1,nocc,
 
770
     &   g_vectors,1,nocc,1,naos,
 
771
     &   g_work   ,1,naos,1,naos)
 
772
      call ga_sync()
 
773
 
 
774
      if (debug) write (6,*) '4'
 
775
 
 
776
c     add the two terms together and transpose the density matrix
 
777
             if (ga_nodeid().eq.0)
 
778
     &       write(*,*) '---- CC1^t-------- START'
 
779
             call ga_print(g_ptmp)
 
780
            if (ga_nodeid().eq.0)
 
781
     &       write(*,*) '---- CC1^t--------- END'           
 
782
             if (ga_nodeid().eq.0)
 
783
     &       write(*,*) '---- C1C^t-------- START'
 
784
             call ga_print(g_work)
 
785
            if (ga_nodeid().eq.0)
 
786
     &       write(*,*) '---- C1C^t--------- END'           
 
787
 
 
788
      call ga_add(1d0,g_ptmp,1d0,g_work,g_work)
 
789
 
 
790
c             if (ga_nodeid().eq.0)
 
791
c     &       write(*,*) '---- g_pmats-0-------- START'
 
792
c             call ga_print(g_work)
 
793
c            if (ga_nodeid().eq.0)
 
794
c     &       write(*,*) '---- g_pmats-0-------- END'           
 
795
 
 
796
      call ga_sync()
 
797
      call ga_transpose(g_work, g_ptmp)
 
798
      call ga_sync()
 
799
 
 
800
c     calculate symmetrized and antisymmetrized part (+ component)
 
801
c     if requested on input:
 
802
 
 
803
      if (lantisym) then
 
804
        call ga_transpose(g_ptmp,g_work)
 
805
        call ga_sync()
 
806
        call ga_add(half,g_ptmp,half,g_work,g_pmats(1))
 
807
        call ga_sync()
 
808
        call ga_add(half,g_ptmp,-half,g_work,g_pmata(1))
 
809
        call ga_sync()
 
810
      else
 
811
        call ga_copy(g_ptmp, g_pmats(1))
 
812
        call ga_sync()
 
813
      endif
 
814
             if (ga_nodeid().eq.0)
 
815
     &       write(*,*) '---- g_pmats-------- START'
 
816
             call ga_print(g_pmats(1))
 
817
            if (ga_nodeid().eq.0)
 
818
     &       write(*,*) '---- g_pmats-------- END'      
 
819
             if (ga_nodeid().eq.0)
 
820
     &       write(*,*) '---- g_pmata-------- START'
 
821
             call ga_print(g_pmata(1))
 
822
            if (ga_nodeid().eq.0)
 
823
     &       write(*,*) '---- g_pmata-------- END'        
 
824
      if (ga_nodeid().eq.0)
 
825
     & write(*,*) 'FA-BEF-get_PertDens'
 
826
      if (ga_nodeid().eq.0) then
 
827
       write(*,1000) ip,im,lstatic,imag,
 
828
     &               haveocc,lantisym,ncomp,naos,nmo,nocc,nvir
 
829
 1000  format('(ip,im,lstatic,imag,haveocc,lanti)=(',
 
830
     &       i2,',',i2,',',L1,',',L1,',',L1,',',L1,') ',
 
831
     &       ' (ncomp,naos,nmo,nocc,nvir)=(',
 
832
     &       i4,',',i4,',',i4,',',i4,',',i4,')')
 
833
      endif
 
834
 
 
835
      call get_PertDens(
 
836
     &               g_pmats(1),!out:      symmetrized dens
 
837
     &               g_pmata(1),!out: anti symmetrized dens
 
838
     &               g_amat,    ! in: the u mat
 
839
     &               g_vectors, ! in: MO vect
 
840
     &               ncomp,     ! in: nr. components g_amat
 
841
     &               ip,im,     ! in: indices of u vect
 
842
     &               imag,      ! in: = T -> imag
 
843
     &               haveocc,   ! in: logical flag
 
844
     &               lantisym,  ! in: logical flag
 
845
     &               naos,nmo,  ! in: nr. AOs,MOs
 
846
     &               nocc,nvir, ! in: nr. (occ,virt) MOs
 
847
     &               debug)     ! in: =.true. -> show debug printouts
 
848
 
 
849
             if (ga_nodeid().eq.0)
 
850
     &       write(*,*) '---- g_pmats-nw-------- START'
 
851
             call ga_print(g_pmats(1))
 
852
            if (ga_nodeid().eq.0)
 
853
     &       write(*,*) '---- g_pmats-nw-------- END'      
 
854
             if (ga_nodeid().eq.0)
 
855
     &       write(*,*) '---- g_pmata-nw-------- START'
 
856
             call ga_print(g_pmata(1))
 
857
            if (ga_nodeid().eq.0)
 
858
     &       write(*,*) '---- g_pmata-nw-------- END'        
 
859
      if (ga_nodeid().eq.0)
 
860
     & write(*,*) 'FA-AFT-get_PertDens'
 
861
c ====> WARNING-- START
 
862
c -- Commenting lines below to enter 2nd part of calc
 
863
      if (ga_nodeid().eq.0) write(*,*) 'SKIP lstatic-if'
 
864
c ====> WARNING-- END
 
865
c      if (lstatic .or. ncomp.eq.1) then
 
866
c       skip calculation of component 2 of the density matrix
 
867
c        if (ga_nodeid().eq.0)
 
868
c     &   write(*,*) 'FA-Skipping calc of 2nd component'
 
869
 
 
870
c        goto 7000
 
871
c      endif  
 
872
      
 
873
        if (ga_nodeid().eq.0)
 
874
     &   write(*,*) 'FA-Doing calc of 2nd component'     
 
875
 
 
876
c     -----------------------------------------
 
877
c     Next step: assemble P(-). Same as before,
 
878
c     but +/- interchanged in amat:
 
879
c     -----------------------------------------
 
880
 
 
881
      call ga_zero(g_ptmp)
 
882
      call ga_zero(g_work)
 
883
      call ga_zero(g_eig1)
 
884
      call ga_sync()
 
885
 
 
886
c     ----------------------------
 
887
c     First step:  C n C(+,dagger)
 
888
c     ----------------------------
 
889
 
 
890
c     calculate C(+,dagger)
 
891
      two = 2d0
 
892
      if (imag) two = -2d0
 
893
      if (.not.haveocc) then
 
894
      call ga_matmul_patch('n','n', two,0d0,
 
895
     &   g_vectors ,1,naos,nocc+1,nmo,
 
896
     &   g_amat(ip),1,nvir,1,nocc,
 
897
     &   g_eig1    ,1,naos,1,nocc)
 
898
      else
 
899
      call ga_matmul_patch('n','n', two,0d0,
 
900
     &   g_vectors ,1,naos,1,nmo,
 
901
     &   g_amat(ip),1,nmo,1,nocc,
 
902
     &   g_eig1    ,1,naos,1,nocc)
 
903
      endif
 
904
      call ga_sync()
 
905
 
 
906
      if (debug) write (6,*) '5'
 
907
 
 
908
c     calculate C(0)C(+,dagger), store in g_ptmp
 
909
      call ga_matmul_patch('n','t', 1d0,0d0,
 
910
     &   g_vectors,1,naos,1,nocc,
 
911
     &   g_eig1,1,nocc,1,naos,
 
912
     &   g_ptmp,1,naos,1,naos)
 
913
      call ga_sync()
 
914
 
 
915
      if (debug) write (6,*) '6'
 
916
      
 
917
c     -----------------------------
 
918
c     Second step: C(-) n C(dagger)
 
919
c     -----------------------------
 
920
 
 
921
c     calculate C(-)
 
922
      two = 2d0
 
923
      if (.not.haveocc) then
 
924
      call ga_matmul_patch('n','n', two,0d0,
 
925
     &   g_vectors,1,naos,nocc+1,nmo,
 
926
     &   g_amat(im),1,nvir,1,nocc,
 
927
     &   g_eig1,1,naos,1,nocc)
 
928
      else
 
929
      call ga_matmul_patch('n','n', two,0d0,
 
930
     &   g_vectors,1,naos,1,nmo,
 
931
     &   g_amat(im),1,nmo,1,nocc,
 
932
     &   g_eig1,1,naos,1,nocc)
 
933
      endif
 
934
      call ga_sync()
 
935
 
 
936
      if (debug) write (6,*) '7'
 
937
 
 
938
c     calculate C(-)C(0,dagger), store in g_work
 
939
      call ga_matmul_patch('n','t', 1d0,0d0,
 
940
     &   g_eig1,1,naos,1,nocc,
 
941
     &   g_vectors,1,nocc,1,naos,
 
942
     &   g_work,1,naos,1,naos)
 
943
      call ga_sync()
 
944
 
 
945
      if (debug) write (6,*) '8'
 
946
 
 
947
c     add the two terms together and transpose
 
948
      call ga_add(1d0,g_ptmp,1d0,g_work,g_work)
 
949
      call ga_sync()
 
950
      call ga_transpose(g_work, g_ptmp)
 
951
      call ga_sync()
 
952
 
 
953
c     calculate symmetrized and antisymmetrized part (- component)
 
954
 
 
955
      if (lantisym) then
 
956
        call ga_transpose(g_ptmp,g_work)
 
957
        call ga_sync()
 
958
        call ga_add(half,g_ptmp, half,g_work,g_pmats(2))
 
959
        call ga_sync()
 
960
        call ga_add(half,g_ptmp,-half,g_work,g_pmata(2))
 
961
        call ga_sync()
 
962
      else
 
963
        call ga_copy(g_ptmp, g_pmats(2))
 
964
        call ga_sync()
 
965
      endif
 
966
             if (ga_nodeid().eq.0)
 
967
     &       write(*,*) '---- g_pmats-2-------- START'
 
968
             call ga_print(g_pmats(2))
 
969
            if (ga_nodeid().eq.0)
 
970
     &       write(*,*) '---- g_pmats-2-------- END'      
 
971
             if (ga_nodeid().eq.0)
 
972
     &       write(*,*) '---- g_pmata-2-------- START'
 
973
             call ga_print(g_pmata(2))
 
974
            if (ga_nodeid().eq.0)
 
975
     &       write(*,*) '---- g_pmata-2-------- END' 
 
976
      if (ga_nodeid().eq.0)
 
977
     & write(*,*) 'FA-BEF-get_PertDens-2'
 
978
      if (ga_nodeid().eq.0) then
 
979
       write(*,1001) ip,im,lstatic,imag,
 
980
     &               haveocc,lantisym,ncomp,naos,nmo,nocc,nvir
 
981
 1001  format('(ip,im,lstatic,imag,haveocc,lanti)=(',
 
982
     &       i2,',',i2,',',L1,',',L1,',',L1,',',L1,') ',
 
983
     &       ' (ncomp,naos,nmo,nocc,nvir)=(',
 
984
     &       i4,',',i4,',',i4,',',i4,',',i4,')')
 
985
      endif
 
986
 
 
987
      call get_PertDens(
 
988
     &               g_pmats(2),!out:      symmetrized dens
 
989
     &               g_pmata(2),!out: anti symmetrized dens
 
990
     &               g_amat,    ! in: the u mat
 
991
     &               g_vectors, ! in: MO vect
 
992
     &               ncomp,     ! in: nr. components g_amat
 
993
     &               im,ip,     ! in: indices of u vect
 
994
     &               imag,      ! in: = T -> imag
 
995
     &               haveocc,   ! in: logical flag
 
996
     &               lantisym,  ! in: logical flag
 
997
     &               naos,nmo,  ! in: nr. AOs,MOs
 
998
     &               nocc,nvir, ! in: nr. (occ,virt) MOs
 
999
     &               debug)     ! in: =.true. -> show debug printouts
 
1000
 
 
1001
             if (ga_nodeid().eq.0)
 
1002
     &       write(*,*) '---- g_pmats-2-nw-------- START'
 
1003
             call ga_print(g_pmats(2))
 
1004
            if (ga_nodeid().eq.0)
 
1005
     &       write(*,*) '---- g_pmats-2-nw-------- END'      
 
1006
             if (ga_nodeid().eq.0)
 
1007
     &       write(*,*) '---- g_pmata-2-nw-------- START'
 
1008
             call ga_print(g_pmata(2))
 
1009
            if (ga_nodeid().eq.0)
 
1010
     &       write(*,*) '---- g_pmata-2-nw-------- END'        
 
1011
      if (ga_nodeid().eq.0)
 
1012
     & write(*,*) 'FA-AFT-get_PertDens'
 
1013
 
 
1014
      stop
 
1015
      
 
1016
c     ---------------------------------------------
 
1017
c     deallocate temporary arrays, sync, and return
 
1018
c     ---------------------------------------------
 
1019
 
 
1020
c     jump here from above in case of static calculation
 
1021
 7000 continue
 
1022
 
 
1023
      if (ga_nodeid().eq.0)
 
1024
     & write(*,*) 'FA-AFT-get_PertDens-did-not-enter-2'
 
1025
 
 
1026
      stop      
 
1027
 
 
1028
      
 
1029
      if (.not.ga_destroy(g_ptmp))
 
1030
     &   call 
 
1031
     &   errquit('CalcPerturbedTDPmat: ga_destroy failed g_ptmp',
 
1032
     &   0,GA_ERR)
 
1033
      
 
1034
      if (.not.ga_destroy(g_work))
 
1035
     &   call 
 
1036
     &   errquit('CalcPerturbedTDPmat: ga_destroy failed g_work',
 
1037
     &   0,GA_ERR)
 
1038
      
 
1039
      if (.not.ga_destroy(g_eig1))
 
1040
     &   call 
 
1041
     &   errquit('CalcPerturbedTDPmat: ga_destroy failed g_eig1',
 
1042
     &   0,GA_ERR)
 
1043
 
 
1044
      call ga_sync()
 
1045
 
 
1046
c     ==================================================================
 
1047
      
 
1048
      return
 
1049
      end
 
1050
 
 
1051
c ++++++++++++++++++ Added for debugging ++++++++++++++++ END