← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/property/hnd_gshift_zora.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/hnd_gshift_zora.F
82
82
      integer debug_gshift
83
83
      integer shift,disp
84
84
      logical skip_cphf_ev_gshift
85
 
 
86
 
c      switch_skip_cphf=.true. ! For skipping cphf or cpks
 
85
c
 
86
      debug_gshift=0 ! =1 for debugging print outs of matrices
 
87
c     switch_skip_cphf=.true. ! For skipping cphf or cpks
87
88
      switch_skip_cphf=.false. ! For NOT skipping cphf or cpks
88
89
      switch_gshift_analysis=.false. ! For default mode gshift tensor
89
 
c      switch_gshift_analysis=.true.  ! For analysis of gshift tensor
90
 
 
91
 
      debug_gshift=0 ! =1 for debugging print outs of matrices
 
90
c     switch_gshift_analysis=.true.  ! For analysis of gshift tensor
 
91
c
92
92
      npol=2 ! gshift works ONLY for npol=2
93
93
      if (ga_nodeid().eq.0) then
94
94
        write(LuOut,*)
108
108
      do ifld = 1,nind_jk
109
109
        jfac(ifld) =  0.0d0       ! used in shell_fock_build()
110
110
        kfac(ifld) = -1.0d0*xfac  ! used in shell_fock_build()
111
 
c       if (ga_nodeid().eq.0) then
112
 
c       write(*,145) ifld,jfac(ifld),kfac(ifld)
113
 
c 145   format('(j,k)(',i3,')=(',f15.8,',',f15.8,')')
114
 
c       endif
 
111
        if (ga_nodeid().eq.0.and.debug_gshift.eq.1) then
 
112
          write(*,145) ifld,jfac(ifld),kfac(ifld)
 
113
 145   format('(j,k)(',i3,')=(',f15.8,',',f15.8,')')
 
114
        endif
115
115
      enddo
116
116
c     Integral initialization
117
117
      call int_init(rtdb,1,basis)
259
259
     &        '-> closed shell system not allowed!')
260
260
       stop
261
261
      endif 
262
 
c     if (ga_nodeid().eq.0) 
263
 
c    &  write(*,*) 'coeffpol=',coeffpol
 
262
      if (ga_nodeid().eq.0) 
 
263
     &  write(*,*) 'coeffpol=',coeffpol
264
264
 
265
265
c ------ Store nopen in rtdb so that CPHF routine is happy ---- START
266
266
          if (.not. rtdb_put(rtdb, 'scf:nopen', 
286
286
 
287
287
c ----- TEST: for testing Coulomb and Exchange contrib --- START
288
288
      if (switch_gshift_analysis) then
289
 
       if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_Coul))
 
289
       if(.not.ga_create(MT_DBL,ntot,3,'rhs_Coul',-1,-1,g_rhs_Coul))
290
290
     &   call errquit('hnd_gshift: ga_create failed g_rhsJ',0,GA_ERR)
291
291
       call ga_zero(g_rhs_Coul)
292
 
       if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_Exch))
 
292
       if(.not.ga_create(MT_DBL,ntot,3,'rhs_Exch',-1,-1,g_rhs_Exch))
293
293
     &   call errquit('hnd_gshift: ga_create failed g_rhsK',0,GA_ERR)
294
294
       call ga_zero(g_rhs_Exch)
295
 
       if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_noJK))
 
295
       if(.not.ga_create(MT_DBL,ntot,3,'rhs_noJK',-1,-1,g_rhs_noJK))
296
296
     &   call errquit('hnd_gshift: ga_create failed g_rhsnJK',0,GA_ERR)
297
297
       call ga_zero(g_rhs_noJK)
298
 
       if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_eSji))
 
298
       if(.not.ga_create(MT_DBL,ntot,3,'rhs_eSji',-1,-1,g_rhs_eSji))
299
299
     &   call errquit('hnd_gshift: ga_create failed g_rhseSji',0,GA_ERR)
300
300
       call ga_zero(g_rhs_eSji)
301
 
       if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_1e))
 
301
       if(.not.ga_create(MT_DBL,ntot,3,'rhs_1e',-1,-1,g_rhs_1e))
302
302
     &   call errquit('hnd_gshift: ga_create failed g_rhs1e',0,GA_ERR)
303
303
       call ga_zero(g_rhs_1e)
304
304
      endif
492
492
      call ga_zero(g_fock)
493
493
c -- TEST: create g_fock_Coul,g_fock_Exch ----------- START
494
494
      if (switch_gshift_analysis) then
495
 
       if (.not.nga_create(MT_DBL,3,ahi,'Fock matrix',
 
495
       if (.not.nga_create(MT_DBL,3,ahi,'Fock matrix g_fock_Coul',
496
496
     &    alo,g_fock_Coul)) call 
497
497
     &    errquit('hnd_gshift: nga_create failed g_fockJ',0,GA_ERR)
498
498
       call ga_zero(g_fock_Coul)
499
 
       if (.not.nga_create(MT_DBL,3,ahi,'Fock matrix',
 
499
       if (.not.nga_create(MT_DBL,3,ahi,'Fock matrix g_fock_Exch',
500
500
     &    alo,g_fock_Exch)) call 
501
501
     &    errquit('hnd_gshift: nga_create failed g_fockX',0,GA_ERR)
502
502
       call ga_zero(g_fock_Exch)
594
594
        ahi(3) = shift+3
595
595
        if      (npol.eq.1) then
596
596
         call nga_scale_patch(g_rhs_noJK,blo,bhi,-4.0d0)
597
 
         call nga_scale_patch(g_rhs_1e,blo,bhi,-4.0d0)
 
597
         call nga_scale_patch(g_rhs_1e  ,blo,bhi,-4.0d0)
598
598
         call nga_scale_patch(g_rhs_eSji,blo,bhi,-4.0d0)
599
599
        else if (npol.eq.2) then
600
600
         call nga_scale_patch(g_rhs_noJK,blo,bhi,-2.0d0)
601
 
         call nga_scale_patch(g_rhs_1e,blo,bhi,-2.0d0)
 
601
         call nga_scale_patch(g_rhs_1e  ,blo,bhi,-2.0d0)
602
602
         call nga_scale_patch(g_rhs_eSji,blo,bhi,-2.0d0)
603
603
        endif
604
604
       enddo ! end-loop-ispin
612
612
       call add_fock1(g_rhs_Coul, ! out: accumulated rhs expression
613
613
     &                g_fock_Coul,! in 
614
614
     &                nmo,npol,nocc,nvirt) 
 
615
        if (.not.ga_destroy(g_fock_Coul)) call 
 
616
     &  errquit('hnd_gshift_zora: ga_destroy failed g_fock_Coul',
 
617
     &           0,GA_ERR)
615
618
       call add_fock1(g_rhs_Exch, ! out: accumulated rhs expression
616
619
     &                g_fock_Exch,! in 
617
620
     &                nmo,npol,nocc,nvirt) 
 
621
        if (.not.ga_destroy(g_fock_Exch)) call 
 
622
     &  errquit('hnd_gshift_zora: ga_destroy failed g_fock_Exch',
 
623
     &           0,GA_ERR)
618
624
      endif
619
625
      if (debug_gshift.eq.1) then
620
626
       if (ga_nodeid().eq.0) then
641
647
     &                    g_rhs_1e,
642
648
     &                    g_rhs_eSji,
643
649
     &                    dbl_mb(k_eval), ! IN: energy eigenvalues
644
 
     &                    nocc,nvirt,nmo,npol)
 
650
     &                    nocc,
 
651
     &                    nvirt,
 
652
     &                    nbf,  ! FA-04-28-12: replacing nmo by nbf
 
653
     &                    npol)
645
654
       else
646
655
         if (ga_nodeid().eq.0)
647
656
     &     write(*,*) 'WARNING : Using skip_cphf ...'
648
657
           call skip_cphf(g_rhs,          ! IN/OUT
649
658
     &                    dbl_mb(k_eval), ! IN: energy eigenvalues
650
 
     &                    nocc,nvirt,nmo,npol)
 
659
     &                    nocc,
 
660
     &                    nvirt,
 
661
     &                    nbf,  ! FA-04-28-12: replacing nmo by nbf
 
662
     &                    npol)
651
663
       endif
652
664
      endif ! END-skip-cphf
653
665
c ----> SKIPPING CPHF ------------------ END
671
683
     &  skip_cphf_ev_gshift = .false.       
672
684
      if (.not.(switch_skip_cphf) .and.
673
685
     &    .not.(skip_cphf_ev_gshift)) then ! START-check-skip-cphf
674
 
       if (ga_nodeid().eq.0)
675
 
     &  write(*,*) 'COMPUTE cphf g-shift data ...'
 
686
c       if (ga_nodeid().eq.0)
 
687
c     &  write(*,*) 'COMPUTE cphf g-shift data ...'
676
688
c     Write ga_rhs to disk 
677
689
       call util_file_name('cphf_rhs',.true.,.true.,cphf_rhs)
678
690
       call util_file_name('cphf_sol',.true.,.true.,cphf_sol)
816
828
     &          lbl_nlmogshift, ! in: for g-shift NLMO analysis
817
829
     &          lbl_nlmoshield, ! in: for shield  NLMO analysis
818
830
     &          rtdb) 
 
831
        if (.not.ga_destroy(g_rhs_Coul)) call 
 
832
     &  errquit('hnd_gshift_zora: ga_destroy failed g_rhs_Coul',
 
833
     &           0,GA_ERR)
 
834
        if (.not.ga_destroy(g_rhs_Exch)) call 
 
835
     &  errquit('hnd_gshift_zora: ga_destroy failed g_rhs_Exch',
 
836
     &           0,GA_ERR)
 
837
        if (.not.ga_destroy(g_rhs_noJK)) call 
 
838
     &  errquit('hnd_gshift_zora: ga_destroy failed g_rhs_noJK',
 
839
     &           0,GA_ERR)
 
840
        if (.not.ga_destroy(g_rhs_1e)) call 
 
841
     &  errquit('hnd_gshift_zora: ga_destroy failed g_rhs_1e',
 
842
     &           0,GA_ERR)
 
843
        if (.not.ga_destroy(g_rhs_eSji)) call 
 
844
     &  errquit('hnd_gshift_zora: ga_destroy failed g_rhs_eSji',
 
845
     &           0,GA_ERR)
819
846
       else
820
847
c      call get_P10( g_d1, ! out: Perturbed density matrix
821
848
c     &             g_rhs, ! in: accumulated rhs expression
1031
1058
       if (.not.ga_destroy(g_rhs0)) call 
1032
1059
     &    errquit('hnd_gshift_zora: ga_destroy failed g_rhs0',
1033
1060
     &            0,GA_ERR)
1034
 
       do ispin=1,npol
 
1061
       do ispin=1,ndens
1035
1062
        if (.not.ga_destroy(g_dens(ispin))) call 
1036
1063
     &    errquit('hnd_gshift_zora: ga_destroy failed g_dens',
1037
1064
     &    0,GA_ERR)
1216
1243
     &                     eval,  ! IN: energy eigenvalues
1217
1244
     &                     nocc,  ! IN: nr. occupied MOs
1218
1245
     &                     nvirt, ! IN: nr. virtual  MOs
1219
 
     &                     nmo,   ! IN: nr. MOs
 
1246
     &                     nbf,   ! IN: nr. basis functions
1220
1247
     &                     npol)  ! IN: nr. polarizations  
1221
1248
c     Purpose : Avoid usage of cphf() routine in
1222
1249
c               the case of BP functionals.
1226
1253
c               done by routine get_P10_1()
1227
1254
c     Author : Fredy W. Aquino
1228
1255
c     Date   : 02-22-11
 
1256
c     Update-04-28-12: Replacing nmo by nbf
1229
1257
c -------- Calculation g_rhs/(ej-ei) ----- START
1230
1258
      implicit none
1231
1259
#include "errquit.fh"
1236
1264
      integer g_rhs
1237
1265
      integer disp,disp1,iocc,ivirt,ind
1238
1266
      integer ispin,alo(3),ahi(3) 
1239
 
      integer nmo,npol,nocc(2),nvirt(2)
1240
 
      double precision eval(nmo*npol),tosclji,coeff
 
1267
      integer nbf,npol,nocc(2),nvirt(2)
 
1268
      double precision eval(nbf*npol),tosclji,coeff
1241
1269
 
1242
1270
      if (npol.eq.1) then
1243
1271
        coeff=0.25d0
1255
1283
 
1256
1284
       if (ga_nodeid().eq.0) then
1257
1285
        write(*,97) npol,nocc(1),nocc(2),
1258
 
     &              nvirt(1),nvirt(2),nmo    
1259
 
 97     format('(npol,nocc,nvirt,nmo)=[',i5,',(',
 
1286
     &              nvirt(1),nvirt(2),nbf    
 
1287
 97     format('(npol,nocc,nvirt,nbf)=[',i5,',(',
1260
1288
     &          i8,',',i8,'),(',i8,',',i8,'),',
1261
1289
     &          i8,']')               
1262
1290
       endif
1263
1291
       do ispin=1,npol 
1264
 
         disp=nmo*(ispin-1)
 
1292
         disp=nbf*(ispin-1) 
1265
1293
         disp1=nocc(1)*nvirt(1)*(ispin-1)
1266
1294
         do ivirt = 1, nvirt(ispin)
1267
1295
          do iocc = 1, nocc(ispin)
1302
1330
     &                     eval,  ! IN: energy eigenvalues
1303
1331
     &                     nocc,  ! IN: nr. occupied MOs
1304
1332
     &                     nvirt, ! IN: nr. virtual  MOs
1305
 
     &                     nmo,   ! IN: nr. MOs
 
1333
     &                     nbf,   ! IN: nr. basis functions
1306
1334
     &                     npol)  ! IN: nr. polarizations  
1307
1335
c     Purpose : Avoid usage of cphf() routine in
1308
1336
c               the case of BP functionals.
1323
1351
     &        g_rhs_1e,g_rhs_eSji
1324
1352
      integer disp,disp1,iocc,ivirt,ind
1325
1353
      integer ispin,alo(3),ahi(3) 
1326
 
      integer nmo,npol,nocc(2),nvirt(2)
1327
 
      double precision eval(nmo*npol),tosclji,coeff
 
1354
      integer nbf,npol,nocc(2),nvirt(2)
 
1355
      double precision eval(nbf*npol),tosclji,coeff
1328
1356
      logical debug_skip
1329
1357
 
1330
1358
      debug_skip= .false. ! =1 for debugging 
1345
1373
 
1346
1374
       if (ga_nodeid().eq.0) then
1347
1375
        write(*,97) npol,nocc(1),nocc(2),
1348
 
     &              nvirt(1),nvirt(2),nmo    
1349
 
 97     format('(npol,nocc,nvirt,nmo)=[',i5,',(',
 
1376
     &              nvirt(1),nvirt(2),nbf    
 
1377
 97     format('(npol,nocc,nvirt,nbf)=[',i5,',(',
1350
1378
     &          i8,',',i8,'),(',i8,',',i8,'),',
1351
1379
     &          i8,']')               
1352
1380
       endif
1353
1381
       do ispin=1,npol 
1354
 
         disp=nmo*(ispin-1)
 
1382
         disp=nbf*(ispin-1)
1355
1383
         disp1=nocc(1)*nvirt(1)*(ispin-1)
1356
1384
         do ivirt = 1, nvirt(ispin)
1357
1385
          do iocc = 1, nocc(ispin)
1454
1482
      do ispin=1,npol
1455
1483
        ntot=ntot+nocc(ispin)*nocc(ispin)
1456
1484
      enddo
1457
 
      if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs0))
 
1485
      if(.not.ga_create(MT_DBL,ntot,3,'rhs0',-1,-1,g_rhs0))
1458
1486
     &   call errquit('get_prelim_fock: ga_create failed g_rhs0',
1459
1487
     &                 0,GA_ERR)
1460
1488
      call ga_zero(g_rhs0)
1494
1522
     &                    pos(1),natoms,oskel)    
1495
1523
 
1496
1524
       blo(1) = 1
1497
 
       bhi(1) = nmo
 
1525
       bhi(1) = nbf ! nmo-fixing4lineardependency
1498
1526
       blo(2) = 1
1499
 
       bhi(2) = nmo
 
1527
       bhi(2) = nbf ! nmo-fixing4lineardependency
1500
1528
       blo(3) = 1
1501
1529
       bhi(3) = 3
1502
1530
      do ispin=1,npol  
1503
1531
       disp=3*(ispin-1) 
1504
1532
       alo(1) = 1
1505
 
       ahi(1) = nmo
 
1533
       ahi(1) = nbf ! nmo-fixing4lineardependency
1506
1534
       alo(2) = 1
1507
 
       ahi(2) = nmo
 
1535
       ahi(2) = nbf ! nmo-fixing4lineardependency
1508
1536
       alo(3) = disp+1
1509
1537
       ahi(3) = disp+3   
1510
1538
       call nga_copy_patch('n',g_s10_1,blo,bhi,
1521
1549
 
1522
1550
      call giao_aotomo(g_s10,vectors,nocc,nvirt,npol,3,nbf)
1523
1551
      if (debug_prelim.eq.1) then
1524
 
       if (ga_nodeid().eq.0) 
1525
 
     &  write(*,*) ' ----------- g_d2------------ START'
1526
1552
        do ispin=1,npol
 
1553
          if (ga_nodeid().eq.0) 
 
1554
     &     write(*,*) '-------- MO vect(',ispin,')----START'
1527
1555
         call ga_print(vectors(ispin))
 
1556
          if (ga_nodeid().eq.0) 
 
1557
     &     write(*,*) '-------- MO vect(',ispin,')----END'
1528
1558
        enddo
 
1559
       if (ga_nodeid().eq.0) 
 
1560
     &  write(*,*) ' ----gprelim: ---- MO:g_s10 ---- START'
1529
1561
       call ga_print(g_s10)
1530
1562
       if (ga_nodeid().eq.0) 
1531
 
     &  write(*,*) ' ---gprelim: ---- MO:g_s10 ----END'
 
1563
     &  write(*,*) ' ----gprelim: ---- MO:g_s10 ---- END'
1532
1564
      endif
1533
1565
 
1534
1566
      do ispin=1,npol ! ++++++  START-loop-ispin
1557
1589
       ahi(2) = nmo
1558
1590
       alo(3) = disp+1
1559
1591
       ahi(3) = disp+3   
 
1592
       call ga_zero(g_s10_1) ! FA-04-28-12-added
1560
1593
       call nga_copy_patch('n',g_s10,alo,ahi,
1561
1594
     &                       g_s10_1,blo,bhi) 
1562
1595
c     COPY : g_s10 -> g_s10_1 ---------- END    
1578
1611
       ahi(1) = nmo
1579
1612
       alo(3) = 1
1580
1613
       ahi(3) = 3
1581
 
       disp=nmo*(ispin-1)
 
1614
       disp=nbf*(ispin-1)  ! FA-04-28-12-energy-fix
1582
1615
 
1583
1616
       if (debug_prelim.eq.1) then
 
1617
         write(*,*) 'get_prelim_fock: nmo=',nmo,
 
1618
     &              ' nbf=',nbf
1584
1619
         if (ga_nodeid().eq.0) 
1585
1620
     &    write(*,*) ' ---gprelim: ---- BEF-scl:g_s10_1 ----START'
1586
1621
         call ga_print(g_s10_1)
1630
1665
     &                           g_rhs,blo,bhi)
1631
1666
      if (debug_prelim.eq.1) then
1632
1667
        if (ga_nodeid().eq.0) 
1633
 
     &   write(*,*) ' ---gprelim: ---- g_rhs ----START'
 
1668
     &   write(*,*) ' ---gprelim: ---- g_rhs(',ispin,') ----START'
1634
1669
        call ga_print(g_rhs)
1635
1670
        if (ga_nodeid().eq.0) 
1636
 
     &   write(*,*) ' ---gprelim: ---- g_rhs ----END'    
 
1671
     &   write(*,*) ' ---gprelim: ---- g_rhs(',ispin,') ----END'    
1637
1672
      endif
1638
1673
c
1639
1674
c     Construct occ-occ part of the three U matrices
1648
1683
       call nga_scale_patch(g_s10_1,alo,ahi,-0.5d0)
1649
1684
       call nga_copy_patch('n',g_s10_1,alo,ahi,
1650
1685
     &                             g_u,alo,ahi)
 
1686
      if (debug_prelim.eq.1) then
 
1687
        if (ga_nodeid().eq.0) 
 
1688
     &   write(*,*) ' ---gprelim: -- g_u(',ispin,') ----START'
 
1689
        call ga_print(g_u)
 
1690
        if (ga_nodeid().eq.0) 
 
1691
     &   write(*,*) ' ---gprelim: -- g_u(',ispin,') ----END'    
 
1692
      endif
1651
1693
c
1652
1694
c     We also need the occupied-occupied contribution of g_u contributing
1653
1695
c     to the first order density matrix. As this block does not change 
1680
1722
     &                    vectors(ispin),blo,bhi,  
1681
1723
     &                               g_u,alo,ahi,
1682
1724
     &                           g_s10_1,dlo,dhi)  
 
1725
      if (debug_prelim.eq.1) then
 
1726
        if (ga_nodeid().eq.0) 
 
1727
     &   write(*,*) ' ---gprelim: -- g_u-1(',ispin,',',xyz,') ----START'
 
1728
        call ga_print(g_s10_1)
 
1729
        if (ga_nodeid().eq.0) 
 
1730
     &   write(*,*) ' ---gprelim: -- g_u-1(',ispin,',',xyz,') ----END'    
 
1731
      endif
 
1732
        ahi(1) = nocc(ispin)
1683
1733
        ahi(2) = nbf
1684
 
        ahi(1) = nocc(ispin)
1685
1734
        bhi(2) = nocc(ispin)
1686
1735
c     Minus sign as we subtract it from the RHS as we do not include 
1687
1736
c     it in the LHS
1692
1741
     &                        vectors(ispin),blo,bhi,
1693
1742
     &                               g_s10_1,alo,ahi,
1694
1743
     &                                  g_d1,clo,chi)  
 
1744
      if (debug_prelim.eq.1) then
 
1745
        if (ga_nodeid().eq.0) 
 
1746
     &   write(*,*) ' ---gprelim: -- g_d1(',ispin,',',xyz,') ----START'
 
1747
        call ga_print(g_d1)
 
1748
        if (ga_nodeid().eq.0) 
 
1749
     &   write(*,*) ' ---gprelim: -- g_d1(',ispin,',',xyz,') ----END'    
 
1750
      endif
1695
1751
       enddo ! end-loop-xyz
1696
1752
c ------------ back-up g_u --> g_rhs0 ---- START
1697
1753
       alo(1) = 1
1821
1877
       alo(1) = shift+1
1822
1878
       ahi(1) = shift+3   
1823
1879
       alo(2) = 1
1824
 
       ahi(2) = nmo
 
1880
       ahi(2) = nbf ! nmo-fixing4lineardependency
1825
1881
       alo(3) = 1
1826
 
       ahi(3) = nmo
 
1882
       ahi(3) = nbf ! nmo-fixing4lineardependency
1827
1883
       blo(1) = 1
1828
1884
       bhi(1) = 3
1829
1885
       blo(2) = 1
1830
 
       bhi(2) = nmo
 
1886
       bhi(2) = nbf ! nmo-fixing4lineardependency
1831
1887
       blo(3) = 1
1832
 
       bhi(3) = nmo
 
1888
       bhi(3) = nbf ! nmo-fixing4lineardependency
1833
1889
        call nga_copy_patch('n',g_fock,alo,ahi,
1834
1890
     &                           g_fck,blo,bhi) 
1835
1891
c       if (ga_nodeid().eq.0)
1866
1922
       chi(3) = nocc(ispin)
1867
1923
       do xyz = 1,3 ! = x,y,z
1868
1924
         blo(1) = nocc(ispin)+1
 
1925
         bhi(1) = nmo
1869
1926
         blo(2) = 1
1870
 
         bhi(1) = nmo
1871
1927
         bhi(2) = nbf
1872
1928
         alo(3) = xyz
1873
1929
         ahi(3) = xyz
2028
2084
 
2029
2085
      if (.not.(do_zora)) then ! do: HF or DFT (norel calc.)
2030
2086
 
2031
 
        if (ga_nodeid().eq.0)
2032
 
     &   write(*,*) 'enter-NonRel-addH10'
 
2087
c       if (ga_nodeid().eq.0)
 
2088
c    &   write(*,*) 'enter-NonRel-addH10'
2033
2089
 
2034
2090
        call int_giao_1ega(basis,basis,g_s10_1,
2035
2091
     &                     'l10' ,pos(1),natoms,oskel)
2053
2109
     &              1.0d0, ga_Fji, ! update g_s10_1 with ga_Fji
2054
2110
     &              g_s10_1)       ! out
2055
2111
      endif
2056
 
c ++++++++++++ COSMO-added-07-13-11+++++++++++++ START
 
2112
c
2057
2113
c     Get external and cosmo bq contribution
2058
2114
      nbq = 0
2059
2115
      nextbq = 0
2068
2124
        call int_giao_1ega(basis,basis,g_s10_1,'bq10',dbl_mb(k_xyz),
2069
2125
     &                     natoms,oskel)
2070
2126
      end if
2071
 
c ++++++++++++ COSMO-added-07-13-11+++++++++++++ END
 
2127
c
2072
2128
c     Transform H10 to MO and add to g_rhs
2073
2129
c     Copy: g_s10_1 --> g_s10
2074
2130
       blo(1) = 1
2075
 
       bhi(1) = nmo
 
2131
       bhi(1) = nbf ! nmo-fixing4lineardependency
2076
2132
       blo(2) = 1
2077
 
       bhi(2) = nmo
 
2133
       bhi(2) = nbf ! nmo-fixing4lineardependency
2078
2134
       blo(3) = 1
2079
2135
       bhi(3) = 3
2080
2136
      do ispin=1,npol   
2081
2137
       shift=3*(ispin-1)
2082
2138
       alo(1) = 1
2083
 
       ahi(1) = nmo
 
2139
       ahi(1) = nbf ! nmo-fixing4lineardependency
2084
2140
       alo(2) = 1
2085
 
       ahi(2) = nmo
 
2141
       ahi(2) = nbf ! nmo-fixing4lineardependency
2086
2142
       alo(3) = shift+1
2087
2143
       ahi(3) = shift+3   
2088
2144
       call nga_copy_patch('n',g_s10_1,blo,bhi,
2315
2371
       clo(2) = 1
2316
2372
       chi(2) = nbf
2317
2373
       dlo(1) = 1
 
2374
       dhi(1) = nbf
2318
2375
       dlo(2) = 1
2319
 
       dhi(1) = nbf
2320
2376
       dhi(2) = nocc(ispin)
2321
2377
 
2322
2378
c --------- zora scaling of MO vectors(1) ----- START
2327
2383
       call ga_copy(vectors(ispin),vectors_scl(ispin))
2328
2384
       if (do_zora .and. .not.(do_NonRel) .and.
2329
2385
     &     .not.(not_zora_scale)) then
 
2386
 
 
2387
c       if (ga_nodeid().eq.0) write(*,*) 'FA-enter-scaling'
 
2388
c       if (ga_nodeid().eq.0)
 
2389
c     &  write(*,*) '---g_CiFull(',ispin,')------ START'
 
2390
c       call ga_print(g_CiFull(ispin))
 
2391
c       if (ga_nodeid().eq.0)
 
2392
c     &  write(*,*) '---g_CiFull(',ispin,')------ END'
 
2393
 
2330
2394
        call ga_scale_cols(vectors_scl(ispin),g_CiFull(ispin))
2331
2395
       endif
2332
2396
c --------- zora scaling of MO vectors(1) ----- END
2667
2731
     &     errquit('g_p1_ov: nga_create failed g_p1_ov',0,GA_ERR)
2668
2732
       call ga_zero(g_p1_ov)
2669
2733
c ------ For debugging ---- START
2670
 
       if (debug_p10 .eq.1) then
2671
 
        if (.not.nga_create(MT_DBL,3,ahi,'g_p1-b matrix',
2672
 
     &      alo,g_p1)) call 
2673
 
     &     errquit('g_p1: nga_create failed g_p1',0,GA_ERR)
2674
 
        call ga_zero(g_p1)
2675
 
       endif
 
2734
c       if (debug_p10 .eq.1) then
 
2735
c        if (.not.nga_create(MT_DBL,3,ahi,'g_p1-b matrix',
 
2736
c     &      alo,g_p1)) call 
 
2737
c     &     errquit('g_p1: nga_create failed g_p1',0,GA_ERR)
 
2738
c        call ga_zero(g_p1)
 
2739
c       endif
2676
2740
c ------ For debugging ---- END
2677
2741
c +++++++++ store g_u for NMLO analysis +++++++++ START
2678
2742
        if (type_NMR.eq.1 .or. type_NMR.eq.3) then ! store only if Hyperfine calc.
3108
3172
        call ga_print(g_p1_ov)
3109
3173
        if (ga_nodeid().eq.0) 
3110
3174
     &    write(*,*) '------g_p1-2---------- END'
3111
 
       if (ga_nodeid().eq.0) 
3112
 
     &    write(*,*) '------g_p1-tot---------- START'
3113
 
        call ga_add(1.0d0,g_p1_oo,1.0d0,g_p1_ov,g_p1)
3114
 
        call ga_print(g_p1)
3115
 
       if (ga_nodeid().eq.0) 
3116
 
     &    write(*,*) '------g_p1-tot---------- END'
 
3175
c       if (ga_nodeid().eq.0) 
 
3176
c     &    write(*,*) '------g_p1-tot---------- START'
 
3177
c        call ga_add(1.0d0,g_p1_oo,1.0d0,g_p1_ov,g_p1)
 
3178
c        call ga_print(g_p1)
 
3179
c       if (ga_nodeid().eq.0) 
 
3180
c     &    write(*,*) '------g_p1-tot---------- END'
3117
3181
       endif
3118
3182
 
3119
3183
       if (.not.ga_destroy(g_d1_oo)) call 
3689
3753
       ahi(1) = nbf
3690
3754
       ahi(2) = nbf
3691
3755
       ahi(3) = 3
3692
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1_oo)) call 
 
3756
       if (.not.nga_create(MT_DBL,3,ahi,'d1oo matrix',alo,g_d1_oo)) call 
3693
3757
     &     errquit('g_d1_oo: nga_create failed g_d1_oo',0,GA_ERR)
3694
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_p1_oo)) call 
 
3758
       if (.not.nga_create(MT_DBL,3,ahi,'p1oo matrix',alo,g_p1_oo)) call 
3695
3759
     &     errquit('g_p1_oo: nga_create failed g_p1_oo',0,GA_ERR)
3696
3760
       call ga_zero(g_p1_oo)
3697
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1_ov)) call 
 
3761
       if (.not.nga_create(MT_DBL,3,ahi,'d1ov matrix',alo,g_d1_ov)) call 
3698
3762
     &     errquit('g_d1_ov: nga_create failed g_d1_ov',0,GA_ERR)
3699
3763
 
3700
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3764
       if (.not.nga_create(MT_DBL,3,ahi,'d1ov_Coul matrix',
3701
3765
     &     alo,g_d1_ov_Coul)) call 
3702
3766
     &     errquit('g_d1_ovJ: nga_create failed g_d1_ovJ',0,GA_ERR)
3703
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3767
       if (.not.nga_create(MT_DBL,3,ahi,'d1ov_Exch matrix',
3704
3768
     &     alo,g_d1_ov_Exch)) call 
3705
3769
     &     errquit('g_d1_ovK: nga_create failed g_d1_ovK',0,GA_ERR)
3706
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3770
       if (.not.nga_create(MT_DBL,3,ahi,'d1ov_noJK matrix',
3707
3771
     &     alo,g_d1_ov_noJK)) call 
3708
3772
     &     errquit('g_d1_ovnJK: nga_create failed g_d1_ovnJK',
3709
3773
     &     0,GA_ERR)
3710
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3774
       if (.not.nga_create(MT_DBL,3,ahi,'d1ov_1e matrix',
3711
3775
     &     alo,g_d1_ov_1e)) call 
3712
3776
     &     errquit('g_d1_ovnJK: nga_create failed g_d1_ov1e',
3713
3777
     &     0,GA_ERR)
3714
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3778
       if (.not.nga_create(MT_DBL,3,ahi,'d1ov_eSji matrix',
3715
3779
     &     alo,g_d1_ov_eSji)) call 
3716
3780
     &     errquit('g_d1_ovnJK: nga_create failed g_d1_oveSji',
3717
3781
     &     0,GA_ERR)
3718
3782
 
3719
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_p1_ov)) call 
 
3783
       if (.not.nga_create(MT_DBL,3,ahi,'p1ov matrix',alo,g_p1_ov)) call 
3720
3784
     &     errquit('g_p1_ov: nga_create failed g_p1_ov',0,GA_ERR)
3721
3785
       call ga_zero(g_p1_ov)
3722
3786
 
3723
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3787
       if (.not.nga_create(MT_DBL,3,ahi,'p1ov_Coul matrix',
3724
3788
     &     alo,g_p1_ov_Coul)) call 
3725
3789
     &     errquit('g_p1_ovJ: nga_create failed g_p1_ovJ',0,GA_ERR)
3726
3790
       call ga_zero(g_p1_ov_Coul)
3727
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3791
       if (.not.nga_create(MT_DBL,3,ahi,'p1ov_Exch matrix',
3728
3792
     &     alo,g_p1_ov_Exch)) call 
3729
3793
     &     errquit('g_p1_ovK: nga_create failed g_p1_ovK',0,GA_ERR)
3730
3794
       call ga_zero(g_p1_ov_Exch)
3731
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3795
       if (.not.nga_create(MT_DBL,3,ahi,'p1ov_noJK matrix',
3732
3796
     &     alo,g_p1_ov_noJK)) call 
3733
3797
     &     errquit('g_p1_ovnJK: nga_create failed g_p1_ovnJK',
3734
3798
     &     0,GA_ERR)
3735
3799
       call ga_zero(g_p1_ov_noJK)
3736
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3800
       if (.not.nga_create(MT_DBL,3,ahi,'p1ov_1e matrix',
3737
3801
     &     alo,g_p1_ov_1e)) call 
3738
3802
     &     errquit('g_p1_ov1e: nga_create failed g_p1_ov1e',
3739
3803
     &     0,GA_ERR)
3740
3804
       call ga_zero(g_p1_ov_1e)
3741
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
3805
       if (.not.nga_create(MT_DBL,3,ahi,'p1ov_eSji matrix',
3742
3806
     &     alo,g_p1_ov_eSji)) call 
3743
3807
     &     errquit('g_p1_oveSji: nga_create failed g_p1_oveSji',
3744
3808
     &     0,GA_ERR)
3745
3809
       call ga_zero(g_p1_ov_eSji)
3746
3810
c ------ For debugging ---- START
3747
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_p1)) call 
 
3811
       if (.not.nga_create(MT_DBL,3,ahi,'p1 matrix',alo,g_p1)) call 
3748
3812
     &     errquit('g_p1: nga_create failed g_p1',0,GA_ERR)
3749
3813
       call ga_zero(g_p1)
3750
3814
c ------ For debugging ---- END
3893
3957
       call ga_copy(vectors(ispin),vectors_scl(ispin))
3894
3958
       if (do_zora .and. .not.(do_NonRel) .and.
3895
3959
     &     .not.(not_zora_scale)) then
 
3960
 
 
3961
c       if (ga_nodeid().eq.0) write(*,*) 'FA-enter-scaling'
 
3962
c       if (ga_nodeid().eq.0)
 
3963
c     &  write(*,*) '---g_CiFull(',ispin,')------ START'
 
3964
c       call ga_print(g_CiFull(ispin))
 
3965
c       if (ga_nodeid().eq.0)
 
3966
c     &  write(*,*) '---g_CiFull(',ispin,')------ END'
 
3967
 
3896
3968
        call ga_scale_cols(vectors_scl(ispin),g_CiFull(ispin))
3897
3969
       endif
3898
3970
c --------- zora scaling of MO vectors(1) ----- END
4284
4356
       if (ga_nodeid().eq.0) 
4285
4357
     &    write(*,*) '------g_p1-tot---------- END'
4286
4358
       endif
4287
 
 
 
4359
       if (.not.ga_destroy(g_p1)) call 
 
4360
     &    errquit('get_p1: ga_destroy failed g_p1',0,GA_ERR)   
4288
4361
       if (.not.ga_destroy(g_d1_oo)) call 
4289
4362
     &    errquit('get_d1: ga_destroy failed g_d1_oo',0,GA_ERR)   
4290
4363
       if (.not.ga_destroy(g_d1_ov)) call 
4344
4417
     &        dlo(3), dhi(3)
4345
4418
      integer natoms
4346
4419
      double precision pos(3*natoms)
4347
 
      double precision eval(nmo*npol),toscl
 
4420
      double precision eval(nbf*npol),toscl
4348
4421
      integer basis ! basis handle
4349
4422
      integer vectors(npol)
4350
4423
      logical oskel
4351
 
      external int_giao_1ega,giao_aotomo
 
4424
      external int_giao_1ega,
 
4425
     &         giao_aotomo
4352
4426
      integer debug_prelim
4353
4427
      integer ndens
4354
4428
 
4361
4435
        ntot=ntot+nocc(ispin)*nocc(ispin)
4362
4436
      enddo
4363
4437
 
4364
 
      if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs0))
 
4438
      if(.not.ga_create(MT_DBL,ntot,3,'prelim_fock_debug g_rhs0',
 
4439
     &                  -1,-1,g_rhs0))
4365
4440
     &   call errquit('get_prelim_fock: ga_create failed g_rhs0',
4366
4441
     &                 0,GA_ERR)
4367
4442
      call ga_zero(g_rhs0)
4371
4446
      chi(1) = 1  
4372
4447
      chi(2) = -1 
4373
4448
      chi(3) = -1
4374
 
      if (.not.nga_create(MT_DBL,3,clo,'D10 matrix',
 
4449
      if (.not.nga_create(MT_DBL,3,clo,'prelim_fock_debug g_d1',
4375
4450
     &                    chi,g_d1)) 
4376
4451
     &  call errquit('gprelim_fock: nga_create failed g_d1',
4377
4452
     &                0,GA_ERR)
4394
4469
     &                    alo,g_s10)) 
4395
4470
     &  call errquit('gprelim_fock: nga_create failed g_s10',
4396
4471
     &               0,GA_ERR)
 
4472
c ------------ creating ga arrays ---------- END
4397
4473
       call ga_zero(g_s10)
4398
 
 
4399
 
c ------------ creating ga arrays ---------- END
4400
4474
       call int_giao_1ega(basis,basis,g_s10_1,'s10',
4401
 
     &                    pos(1),natoms,oskel)    
 
4475
     &                    pos,natoms,oskel)   
 
4476
      if (debug_prelim.eq.1) then
 
4477
       if (ga_nodeid().eq.0) 
 
4478
     &  write(*,*) ' ---gprelim: ---- g_s10_1 ----START'
 
4479
        call ga_print(g_s10_1)
 
4480
       if (ga_nodeid().eq.0) 
 
4481
     &  write(*,*) ' ---gprelim: ---- g_s10_1 ----END'
 
4482
      endif
4402
4483
 
4403
4484
       blo(1) = 1
4404
 
       bhi(1) = nmo
 
4485
       bhi(1) = nbf ! nmo-fixing4lineardependency
4405
4486
       blo(2) = 1
4406
 
       bhi(2) = nmo
 
4487
       bhi(2) = nbf ! nmo-fixing4lineardependency
4407
4488
       blo(3) = 1
4408
4489
       bhi(3) = 3
4409
4490
      do ispin=1,npol  
4410
4491
       disp=3*(ispin-1) 
4411
4492
       alo(1) = 1
4412
 
       ahi(1) = nmo
 
4493
       ahi(1) = nbf ! nmo-fixing4lineardependency
4413
4494
       alo(2) = 1
4414
 
       ahi(2) = nmo
 
4495
       ahi(2) = nbf ! nmo-fixing4lineardependency
4415
4496
       alo(3) = disp+1
4416
4497
       ahi(3) = disp+3   
4417
4498
       call nga_copy_patch('n',g_s10_1,blo,bhi,
4425
4506
       if (ga_nodeid().eq.0) 
4426
4507
     &  write(*,*) ' ---gprelim: ---- g_s10 ----END'
4427
4508
      endif
4428
 
 
 
4509
c Note.- Output g_s10 is (nmo,nmo) not (nbf,nbf) matrix
4429
4510
      call giao_aotomo(g_s10,vectors,nocc,nvirt,npol,3,nbf)
 
4511
 
4430
4512
      if (debug_prelim.eq.1) then
4431
 
       if (ga_nodeid().eq.0) 
4432
 
     &  write(*,*) ' ----------- g_d2------------ START'
4433
4513
        do ispin=1,npol
 
4514
          if (ga_nodeid().eq.0) 
 
4515
     &     write(*,*) '-------- MO vect(',ispin,')----START'
4434
4516
         call ga_print(vectors(ispin))
 
4517
          if (ga_nodeid().eq.0) 
 
4518
     &     write(*,*) '-------- MO vect(',ispin,')----END'
4435
4519
        enddo
 
4520
       if (ga_nodeid().eq.0) 
 
4521
     &  write(*,*) ' ----gprelim: ---- MO:g_s10 ---- START'
4436
4522
       call ga_print(g_s10)
4437
4523
       if (ga_nodeid().eq.0) 
4438
 
     &  write(*,*) ' ---gprelim: ---- MO:g_s10 ----END'
 
4524
     &  write(*,*) ' ----gprelim: ---- MO:g_s10 ---- END'
4439
4525
      endif
4440
4526
 
4441
4527
      do ispin=1,npol ! ++++++  START-loop-ispin
4463
4549
       alo(2) = 1
4464
4550
       ahi(2) = nmo
4465
4551
       alo(3) = disp+1
4466
 
       ahi(3) = disp+3   
 
4552
       ahi(3) = disp+3  
 
4553
       call ga_zero(g_s10_1) ! FA-04-28-12-added
4467
4554
       call nga_copy_patch('n',g_s10,alo,ahi,
4468
4555
     &                       g_s10_1,blo,bhi) 
4469
4556
c     COPY : g_s10 -> g_s10_1 ---------- END    
4485
4572
       ahi(1) = nmo
4486
4573
       alo(3) = 1
4487
4574
       ahi(3) = 3
4488
 
       disp=nmo*(ispin-1)
 
4575
       disp=nbf*(ispin-1)  ! FA-04-28-12-energy-fix
4489
4576
 
4490
4577
       if (debug_prelim.eq.1) then
4491
4578
         if (ga_nodeid().eq.0) 
4502
4589
 
4503
4590
          if (debug_prelim.eq.1) then
4504
4591
           if (ga_nodeid().eq.0) then
4505
 
            write(*,1) iocc,toscl
4506
 
 1          format('E2scl(',i5,')=',f15.8) 
 
4592
            write(*,1) ispin,iocc,toscl
 
4593
 1          format('E2scl(',i4,',',i5,')=',f15.8) 
4507
4594
           endif
4508
4595
          endif
4509
4596
 
4511
4598
        enddo ! end-loop-iocc
4512
4599
      if (debug_prelim.eq.1) then
4513
4600
        if (ga_nodeid().eq.0) 
4514
 
     &   write(*,*) ' ---gprelim: ---- AFT-scl:g_s10_1 ----START'
 
4601
     &   write(*,*) ' ---gprelim: - AFT-scl:g_s10_1(',ispin,') --START'
4515
4602
        call ga_print(g_s10_1)
4516
4603
        if (ga_nodeid().eq.0) 
4517
 
     &   write(*,*) ' ---gprelim: ---- AFT-scl:g_s10_1 ----END'    
 
4604
     &   write(*,*) ' ---gprelim: - AFT-scl:g_s10_1(',ispin,') --END'    
4518
4605
      endif
4519
4606
 
4520
4607
c     Copy to g_rhs 
4541
4628
c ----- output: g_rhs_eSji --- END
4542
4629
      if (debug_prelim.eq.1) then
4543
4630
        if (ga_nodeid().eq.0) 
4544
 
     &   write(*,*) ' ---gprelim: ---- g_rhs ----START'
 
4631
     &   write(*,*) ' ---gprelim: ---- g_rhs(',ispin,') ----START'
4545
4632
        call ga_print(g_rhs)
4546
4633
        if (ga_nodeid().eq.0) 
4547
 
     &   write(*,*) ' ---gprelim: ---- g_rhs ----END'    
 
4634
     &   write(*,*) ' ---gprelim: ---- g_rhs(',ispin,') ----END'    
4548
4635
      endif
4549
4636
c
4550
4637
c     Construct occ-occ part of the three U matrices
4559
4646
       call nga_scale_patch(g_s10_1,alo,ahi,-0.5d0)
4560
4647
       call nga_copy_patch('n',g_s10_1,alo,ahi,
4561
4648
     &                             g_u,alo,ahi)
 
4649
c         write(*,13) ispin,xyz,nmo,nbf,
 
4650
c     &              alo(1),ahi(1),alo(2),ahi(2),
 
4651
c     &              alo(3),ahi(3)
 
4652
c 13      format('prelim-x: (ispin,xyz,nmo,nbf)=(',
 
4653
c     &          i3,',',i3,',',i3,',',i3,') ',
 
4654
c     &          'alo-ahi=(',i3,',',i3,',',
 
4655
c     &          i3,',',i3,',',i3,',',i3,')')
 
4656
      if (debug_prelim.eq.1) then
 
4657
        if (ga_nodeid().eq.0) 
 
4658
     &   write(*,*) ' ---gprelim: -- g_u(',ispin,') ----START'
 
4659
        call ga_print(g_u)
 
4660
        if (ga_nodeid().eq.0) 
 
4661
     &   write(*,*) ' ---gprelim: -- g_u(',ispin,') ----END'    
 
4662
      endif
4562
4663
c
4563
4664
c     We also need the occupied-occupied contribution of g_u contributing
4564
4665
c     to the first order density matrix. As this block does not change 
4568
4669
        alo(1) = 1
4569
4670
        alo(2) = 1
4570
4671
        blo(1) = 1
 
4672
        bhi(1) = nbf
4571
4673
        blo(2) = 1
4572
 
        bhi(1) = nbf
4573
4674
        clo(2) = 1
 
4675
        chi(2) = nbf
4574
4676
        clo(3) = 1
4575
 
        chi(2) = nbf
4576
4677
        chi(3) = nbf
4577
4678
        dlo(1) = 1
 
4679
        dhi(1) = nbf
4578
4680
        dlo(2) = 1
4579
 
        dhi(1) = nbf
4580
4681
        dhi(2) = nocc(ispin)
4581
4682
c     Create "perturbed density matrix" for closed-closed g_u block
4582
4683
       do xyz = 1,3 ! = x,y,z
4587
4688
        ahi(1) = nmo
4588
4689
        ahi(2) = nocc(ispin)
4589
4690
        bhi(2) = nmo 
 
4691
 
 
4692
c         write(*,3) ispin,xyz,nmo,nbf,
 
4693
c     &              blo(1),bhi(1),blo(2),bhi(2),
 
4694
c     &              blo(3),bhi(3),
 
4695
c     &              alo(1),ahi(1),alo(2),ahi(2),
 
4696
c     &              alo(3),ahi(3),
 
4697
c     &              dlo(1),dhi(1),dlo(2),dhi(2),
 
4698
c     &              dlo(3),dhi(3)
 
4699
c 3      format('prelim-1: (ispin,xyz,nmo,nbf)=(',
 
4700
c     &          i3,',',i3,',',i3,',',i3,') ',
 
4701
c     &          'blo-bhi=(',i3,',',i3,',',
 
4702
c     &          i3,',',i3,',',i3,',',i3,') ',
 
4703
c     &          'alo-ahi=(',i3,',',i3,',',
 
4704
c     &          i3,',',i3,',',i3,',',i3,') ',
 
4705
c     &          'dlo-dhi=(',i3,',',i3,',',
 
4706
c     &          i3,',',i3,',',i3,',',i3,')')
 
4707
 
4590
4708
        call nga_matmul_patch('n','n',1.0d0,0.0d0,
4591
4709
     &                    vectors(ispin),blo,bhi,  
4592
4710
     &                               g_u,alo,ahi,
4593
4711
     &                           g_s10_1,dlo,dhi)  
 
4712
      if (debug_prelim.eq.1) then
 
4713
        if (ga_nodeid().eq.0) 
 
4714
     &   write(*,*) ' ---gprelim: -- g_u-1(',ispin,',',xyz,') ----START'
 
4715
        call ga_print(g_s10_1)
 
4716
        if (ga_nodeid().eq.0) 
 
4717
     &   write(*,*) ' ---gprelim: -- g_u-1(',ispin,',',xyz,') ----END'    
 
4718
      endif
 
4719
        ahi(1) = nocc(ispin)
4594
4720
        ahi(2) = nbf
4595
 
        ahi(1) = nocc(ispin)
4596
4721
        bhi(2) = nocc(ispin)
4597
4722
c     Minus sign as we subtract it from the RHS as we do not include 
4598
4723
c     it in the LHS
4599
4724
        disp=3*(ispin-1)
4600
4725
        clo(1) = disp+xyz
4601
4726
        chi(1) = disp+xyz
 
4727
 
 
4728
c         write(*,4) ispin,xyz,nmo,nbf,
 
4729
c     &              blo(1),bhi(1),blo(2),bhi(2),
 
4730
c     &              blo(3),bhi(3),
 
4731
c     &              alo(1),ahi(1),alo(2),ahi(2),
 
4732
c     &              alo(3),ahi(3),
 
4733
c     &              clo(1),chi(1),clo(2),chi(2),
 
4734
c     &              clo(3),chi(3)
 
4735
c 4      format('prelim-2: (ispin,xyz,nmo,nbf)=(',
 
4736
c     &          i3,',',i3,',',i3,',',i3,') ',
 
4737
c     &          'blo-bhi=(',i3,',',i3,',',
 
4738
c     &          i3,',',i3,',',i3,',',i3,') ',
 
4739
c     &          'alo-ahi=(',i3,',',i3,',',
 
4740
c     &          i3,',',i3,',',i3,',',i3,') ',
 
4741
c     &          'clo-chi=(',i3,',',i3,',',
 
4742
c     &          i3,',',i3,',',i3,',',i3,')')
 
4743
 
 
4744
 
4602
4745
        call nga_matmul_patch('n','t',-1.0d0,0.0d0,
4603
4746
     &                        vectors(ispin),blo,bhi,
4604
4747
     &                               g_s10_1,alo,ahi,
4605
4748
     &                                  g_d1,clo,chi)  
 
4749
      if (debug_prelim.eq.1) then
 
4750
        if (ga_nodeid().eq.0) 
 
4751
     &   write(*,*) ' ---gprelim: -- g_d1(',ispin,',',xyz,') ----START'
 
4752
        call ga_print(g_d1)
 
4753
        if (ga_nodeid().eq.0) 
 
4754
     &   write(*,*) ' ---gprelim: -- g_d1(',ispin,',',xyz,') ----END'    
 
4755
      endif
4606
4756
       enddo ! end-loop-xyz
4607
4757
c ------------ back-up g_u --> g_rhs0 ---- START
4608
4758
       alo(1) = 1
4618
4768
       bhi(2) = 3
4619
4769
       call nga_copy_patch('n',g_u   ,alo,ahi,
4620
4770
     &                         g_rhs0,blo,bhi) ! copy to g_rhs0
4621
 
c ------------ back-up g_u --> g_rhs0 ---- END
4622
 
c ----- Remove scratch ga arrays ------
 
4771
c ---- back-up g_u --> g_rhs0 ---- END
 
4772
c ---- Remove scratch ga arrays ------
4623
4773
       if (.not.ga_destroy(g_u)) call 
4624
4774
     &     errquit('gprelim_fock: ga_destroy failed g_u',0,GA_ERR)
4625
4775
      enddo ! end-loop-ispin
4632
4782
      chi(1) =  1  
4633
4783
      chi(2) = -1 
4634
4784
      chi(3) = -1
4635
 
      if (.not.nga_create(MT_DBL,3,clo,'D10 matrix',
 
4785
      if (.not.nga_create(MT_DBL,3,clo,'prelim_fock_debug g_d2',
4636
4786
     &                    chi,g_d2)) 
4637
4787
     &  call errquit('gprelim_fock: nga_create failed g_d2',
4638
4788
     &                0,GA_ERR)
4724
4874
      logical oskel
4725
4875
      integer NoKinetic
4726
4876
      common /skipKinetic/NoKinetic ! goes to int_giao_1ega()
4727
 
      external int_giao_1ega,giao_aotomo
 
4877
      external int_giao_1ega,
 
4878
     &         giao_aotomo
4728
4879
c     bq charges
4729
4880
      integer nbq,nextbq,ncosbq,
4730
4881
     &        l_xyz,k_xyz ! dummy variables
4752
4903
       call ga_zero(g_s10)
4753
4904
 
4754
4905
      if (.not.(do_zora)) then ! do: HF or DFT (norel calc.)
 
4906
c        call int_giao_1ega(basis,basis,g_s10_1,
 
4907
c     &                     'l10' ,pos(1),natoms,oskel)
4755
4908
        call int_giao_1ega(basis,basis,g_s10_1,
4756
 
     &                     'l10' ,pos(1),natoms,oskel)
 
4909
     &                     'l10' ,pos,natoms,oskel)
4757
4910
        NoKinetic=0 ! =0 DO-kinetic, =1 SKIP-kinetic
 
4911
c        call int_giao_1ega(basis,basis,g_s10_1,
 
4912
c     &                     'tv10',pos(1),natoms,oskel)
4758
4913
        call int_giao_1ega(basis,basis,g_s10_1,
4759
 
     &                     'tv10',pos(1),natoms,oskel)
 
4914
     &                     'tv10',pos,natoms,oskel)
4760
4915
      else                            ! do: zora (relativistic calc.)
4761
4916
c ----------- Create scratch ga-arrays ------- END
4762
4917
        NoKinetic=1 ! =0 DO-kinetic, =1 SKIP-kinetic
 
4918
c        call int_giao_1ega(basis,basis,g_s10_1,
 
4919
c     &                     'tv10',pos(1),natoms,oskel)
4763
4920
        call int_giao_1ega(basis,basis,g_s10_1,
4764
 
     &                     'tv10',pos(1),natoms,oskel)
 
4921
     &                     'tv10',pos,natoms,oskel)
4765
4922
        call ga_add(1.0d0,g_s10_1,
4766
 
     &              1.0d0, ga_Fji, ! update g_s10_1 with ga_Fji
4767
 
     &              g_s10_1)       ! out
 
4923
     &              1.0d0,ga_Fji, ! update g_s10_1 with ga_Fji
 
4924
     &                    g_s10_1)! out
4768
4925
      endif
4769
 
c ++++++++++++ COSMO-added-07-13-11+++++++++++++ START
 
4926
c
4770
4927
c     Get external and cosmo bq contribution
4771
4928
      nbq = 0
4772
4929
      nextbq = 0
4781
4938
        call int_giao_1ega(basis,basis,g_s10_1,'bq10',dbl_mb(k_xyz),
4782
4939
     &                     natoms,oskel)
4783
4940
      end if
4784
 
c ++++++++++++ COSMO-added-07-13-11+++++++++++++ END
 
4941
c
4785
4942
c     Transform H10 to MO and add to g_rhs
4786
4943
c     Copy: g_s10_1 --> g_s10
4787
4944
       blo(1) = 1
4788
 
       bhi(1) = nmo
 
4945
       bhi(1) = nbf ! nmo-fixing4lineardependency
4789
4946
       blo(2) = 1
4790
 
       bhi(2) = nmo
 
4947
       bhi(2) = nbf ! nmo-fixing4lineardependency
4791
4948
       blo(3) = 1
4792
4949
       bhi(3) = 3
4793
4950
      do ispin=1,npol   
4794
4951
       shift=3*(ispin-1)
4795
4952
       alo(1) = 1
4796
 
       ahi(1) = nmo
 
4953
       ahi(1) = nbf ! nmo-fixing4lineardependency
4797
4954
       alo(2) = 1
4798
 
       ahi(2) = nmo
 
4955
       ahi(2) = nbf ! nmo-fixing4lineardependency
4799
4956
       alo(3) = shift+1
4800
4957
       ahi(3) = shift+3   
4801
4958
       call nga_copy_patch('n',g_s10_1,blo,bhi,
4802
4959
     &                           g_s10,alo,ahi) 
4803
4960
      enddo ! end-loop-ispin
4804
4961
      call giao_aotomo(g_s10,vectors,nocc,nvirt,npol,3,nbf)
 
4962
 
4805
4963
      do ispin=1,npol
4806
4964
       shift=3*(ispin-1)
4807
4965
       alo(1) = nocc(ispin)+1
5223
5381
       ahi(1) = nbf
5224
5382
       ahi(2) = nbf
5225
5383
       ahi(3) = 3
5226
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1_oo)) call 
 
5384
       if (.not.nga_create(MT_DBL,3,ahi,'P10_JK_AorB g_d1_oo',
 
5385
     &                     alo,g_d1_oo)) call 
5227
5386
     &     errquit('g_d1_oo: nga_create failed g_d1_oo',0,GA_ERR)
5228
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_p1_oo)) call 
 
5387
       if (.not.nga_create(MT_DBL,3,ahi,'P10_JK_AorB g_p1_oo',
 
5388
     &                     alo,g_p1_oo)) call 
5229
5389
     &     errquit('g_p1_oo: nga_create failed g_p1_oo',0,GA_ERR)
5230
5390
       call ga_zero(g_p1_oo)
5231
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1_ov)) call 
 
5391
       if (.not.nga_create(MT_DBL,3,ahi,'P10_JK_AorB g_d1_ov',
 
5392
     &                     alo,g_d1_ov)) call 
5232
5393
     &     errquit('g_d1_ov: nga_create failed g_d1_ov',0,GA_ERR)
5233
5394
 
5234
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
5395
       if (.not.nga_create(MT_DBL,3,ahi,'P10_JK_AorB g_d1_ov_Coul',
5235
5396
     &     alo,g_d1_ov_Coul)) call 
5236
5397
     &     errquit('g_d1_ovJ: nga_create failed g_d1_ovJ',0,GA_ERR)
5237
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
5398
       if (.not.nga_create(MT_DBL,3,ahi,'P10_JK_AorB g_d1_ovJ',
5238
5399
     &     alo,g_d1_ov_Exch)) call 
5239
5400
     &     errquit('g_d1_ovK: nga_create failed g_d1_ovK',0,GA_ERR)
5240
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
5401
       if (.not.nga_create(MT_DBL,3,ahi,'P10_JK_AorB g_d1_ovnoJK',
5241
5402
     &     alo,g_d1_ov_noJK)) call 
5242
5403
     &     errquit('g_d1_ovnJK: nga_create failed g_d1_ovnJK',
5243
5404
     &     0,GA_ERR)
5272
5433
     &     errquit('g_p1_ov1e: nga_create failed g_p1_ov1e',
5273
5434
     &     0,GA_ERR)
5274
5435
       call ga_zero(g_p1_ov_1e)
5275
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',
 
5436
       if (.not.nga_create(MT_DBL,3,ahi,'P10_JK_AorB g_p1_ov_eSji',
5276
5437
     &     alo,g_p1_ov_eSji)) call 
5277
5438
     &     errquit('g_p1_oveSji: nga_create failed g_p1_oveSji',
5278
5439
     &     0,GA_ERR)
5279
5440
       call ga_zero(g_p1_ov_eSji)
5280
5441
c ------ For debugging ---- START
5281
 
       if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_p1)) call 
5282
 
     &     errquit('g_p1: nga_create failed g_p1',0,GA_ERR)
5283
 
       call ga_zero(g_p1)
 
5442
c       if (.not.nga_create(MT_DBL,3,ahi,
 
5443
c     &          'P10_JK_AorB g_p1',alo,g_p1)) call 
 
5444
c     &     errquit('g_p1: nga_create failed g_p1',0,GA_ERR)
 
5445
c       call ga_zero(g_p1)
5284
5446
c ------ For debugging ---- END
5285
5447
 
5286
5448
c      do ispin=1,npol
5703
5865
        call ga_print(g_p1_ov)
5704
5866
        if (ga_nodeid().eq.0) 
5705
5867
     &    write(*,*) '------g_p1-2---------- END'
5706
 
       if (ga_nodeid().eq.0) 
5707
 
     &    write(*,*) '------g_p1-tot---------- START'
5708
 
        call ga_add(1.0d0,g_p1_oo,1.0d0,g_p1_ov,g_p1)
5709
 
        call ga_print(g_p1)
5710
 
       if (ga_nodeid().eq.0) 
5711
 
     &    write(*,*) '------g_p1-tot---------- END'
 
5868
c       if (ga_nodeid().eq.0) 
 
5869
c     &    write(*,*) '------g_p1-tot---------- START'
 
5870
c        call ga_add(1.0d0,g_p1_oo,1.0d0,g_p1_ov,g_p1)
 
5871
c        call ga_print(g_p1)
 
5872
c       if (ga_nodeid().eq.0) 
 
5873
c     &    write(*,*) '------g_p1-tot---------- END'
5712
5874
       endif
5713
 
 
 
5875
c       if (.not.ga_destroy(g_p1)) call 
 
5876
c     &    errquit('get_p1: ga_destroy failed g_d1_oo',0,GA_ERR)   
5714
5877
       if (.not.ga_destroy(g_d1_oo)) call 
5715
5878
     &    errquit('get_d1: ga_destroy failed g_d1_oo',0,GA_ERR)   
5716
5879
       if (.not.ga_destroy(g_d1_ov)) call 
6205
6368
      return
6206
6369
      end
6207
6370
c -------------- for g-shift NMLO analysis ----------------- END
6208
 
c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6209
 
c $Id: hnd_gshift_zora.F 21439 2011-11-07 23:26:15Z niri $