comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-29 06:35:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

More progress

Browse source
This commit is contained in:
Dave Bernholdt 1998-07-14 21:41:49 +00:00
parent c26cecca5c
commit 75ff08d0b9
4 changed files with 206 additions and 63 deletions
Download patch file Download diff file Expand all files Collapse all files

13
src/rimp2_grad/driver_e.F
View file

@ -3,7 +3,7 @@ C NAME
C RIMP2_Driver_E -- Master routine for RI-MP2 energy evaluation
C
C REVISION
C $Id: driver_e.F,v 1.2 1998-07-13 02:38:07 bernhold Exp $
C $Id: driver_e.F,v 1.3 1998-07-14 21:41:47 bernhold Exp $
C
C SYNOPSIS
Logical Function RIMP2_DRIVER_E(MaxSpin, BraKetSame, FN_Int,
@ -495,17 +495,6 @@ C
Status = GA_Destroy( G_PairE(ISpin+JSpin-1) )
EndDo
EndDo
C
C ****************************************
C * Finished with gradient intermediates *
C ****************************************
C
If ( DoGrad ) then
Do ISpin = 1, TopSpin
If (.NOT. GA_Destroy( G_P2(ISpin) ) ) Call ErrQuit(
$ 'RIMP2_Driver_E: can''t destroy P2 ', 0)
EndDo
EndIf
C
Call GA_Sync
Call PStat_Off(PS_Energy)

28
src/rimp2_grad/driver_g.F
View file

@ -3,11 +3,11 @@ C NAME
C RIMP2_Driver_G -- Master routine for RI-MP2 gradient evaluation
C
C REVISION
C $Id: driver_g.F,v 1.3 1998-07-14 16:18:34 bernhold Exp $
C $Id: driver_g.F,v 1.4 1998-07-14 21:41:48 bernhold Exp $
C
C SYNOPSIS
Logical Function RIMP2_Driver_G(MaxSpin, TopSpin, NFrzO, NAct,
$ NVir, FitBas, FN_Int, FN_Gam, BraKetSame)
$ NVir, FitBas, FN_Int, FN_Gam, BraKetSame, g_P2, Eig, LDEig)
Implicit NONE
C
Integer MaxSpin ![in]
@ -19,6 +19,9 @@ C
Character*(*) FN_Int(MaxSpin, 2, 2) ![in]
Character*(*) FN_Gam(TopSpin) ![in]
Logical BraKetSame ![in]
Integer g_P2(TopSpin) ![in]
Integer LDEig ![in]
Double Precision Eig(LDEig, TopSpin) ![in]
C
C INCLUDE FILES
#include "mafdecls.fh"
@ -38,8 +41,7 @@ C
C
C LOCAL VARIABLES
Integer ISpin, NFit
Integer g_L1(MyMaxSpin), g_L2(MyMaxSpin), g_L3(MyMaxSpin),
$ g_L4(MyMaxSpin), g_L(MyMaxspin)
Integer g_L(MyMaxspin), g_W2(MyMaxSpin)
Logical PrInfo, PrPrgRpt
Character*256 String1
Integer D_Int(MyMaxSpin, 2, 2), D_Gam(MyMaxSpin)
@ -119,10 +121,19 @@ C
$ String1(:Inp_StrLen(String1)), MinChunk,
$ MinChunk, g_L(ISpin) ) ) Call ErrQuit(
$ 'RIMP2_Driver_G: can''t allocate L', ISpin)
C
C Create W2(p,q) (initial contributions come from L)
C
String1 = 'W2 spin ' // SpinItoA(ISpin)
If ( .NOT. GA_Create(MT_Dbl, NFrzO+NAct(ISpin)+NVir(ISpin),
$ NFrzO+NAct(ISpin)+NVir(ISpin),
$ String1(:Inp_StrLen(String1)), MinChunk,
$ MinChunk, g_W2(ISpin) ) ) Call ErrQuit(
$ 'RIMP2_Driver_G: can''t allocate W2', ISpin)
EndDo ! ISpin
C
Call RIMP2_Mk_L(TopSpin, NFrzO, NAct, NVir, NFit, D_Int(1, 1, 1),
$ D_Int(1, 1, 2), D_Gam, g_L)
$ D_Int(1, 1, 2), D_Gam, g_L, g_P2, g_W2, Eig, LDEig)
C
C *******************
C * Clean up memory *
@ -131,6 +142,13 @@ C
Do ISpin = 1, TopSpin
If ( .NOT. GA_Destroy( G_L(ISpin) ) ) Call ErrQuit(
$ 'RIMP2_Driver_G: can''t free L', ISpin)
If ( .NOT. GA_Destroy( G_W2(ISpin) ) ) Call ErrQuit(
$ 'RIMP2_Driver_G: can''t free W2', ISpin)
C
C Created in RIMP2_Driver_E
C
If ( .NOT. GA_Destroy( G_P2(ISpin) ) ) Call ErrQuit(
$ 'RIMP2_Driver_G: can''t free P2', ISpin)
EndDo
C
C ************************

216
src/rimp2_grad/mk_l.F
View file

@ -3,11 +3,11 @@ C NAME
C RIMP2_Mk_L -- Form Lagrangian terms & their contributions elsewhere
C
C REVISION
C $Id: mk_l.F,v 1.1 1998-07-14 16:18:36 bernhold Exp $
C $Id: mk_l.F,v 1.2 1998-07-14 21:41:49 bernhold Exp $
C
C SYNOPSIS
Subroutine RIMP2_Mk_L(TopSpin, NFrzO, NAct, NVir, NFit,
$ D_Int_ai, D_Int_ij, D_Gam, g_L)
$ D_Int_ai, D_Int_ij, D_Gam, g_L, g_P2, g_W2, Eig, LDEig)
Implicit NONE
C
Integer TopSpin ![in]
@ -19,6 +19,10 @@ C
Integer D_Int_ij(TopSpin) ![in]
Integer D_Gam(TopSpin) ![in]
Integer g_L(TopSpin) ![in]
Integer g_P2(TopSpin) ![in]
Integer g_W2(TopSpin) ![in]
Integer LDEig ![in]
Double Precision Eig(LDEig, TopSpin) ![in]
C
C INCLUDE FILES
#include "mafdecls.fh"
@ -39,6 +43,9 @@ C LOCAL VARIABLES
Integer ISpin
Integer g_L1(MaxSpin), g_L2(MaxSpin), g_L3(MaxSpin),g_L4(MaxSpin)
Character*(256) String1
Integer C, A, V, PLo, PHi, QLo, QHi, I, Y, Index, LD
Integer Me
Double Precision Scale
C
C STATEMENT FUNCTIONS
Character*1 SpinItoA
@ -49,6 +56,8 @@ C
If ( TopSpin .gt. MaxSpin) Call ErrQuit(
$ 'RIMP2_Mk_L: fatal program error: TopSpin > MaxSpin',
$ MaxSpin)
C
Me = GA_NodeID()
C
C Allocate matrices for partial results
C
@ -122,68 +131,189 @@ C ********************************
C
Do ISpin = 1, TopSpin
Call GA_Zero( g_L(ISpin) )
C
C Shorthand, so we can actually read the code!
C
C = NFrzO
A = NAct(ISpin)
V = NVir(ISpin)
C
If ( NFrzO .gt. 0) then
C
C L(iy) = -L1(yi)
C
Call GA_Copy('T',
$ g_L1(ISpin), 1, NFrzO, NFrzO+1, NFrzO+NAct(ISpin),
$ g_L(ISpin), NFrzO+1, NFrzO+NAct(ISpin), 1, NFrzO)
Call GA_Scale_Patch(g_L(ISpin), NFrzO+1, NFrzO+NAct(ISpin),
$ 1, NFrzO, -1.0d0 )
Call GA_Copy_Patch('T', g_L1(ISpin), 1, C, C+1, C+A,
$ g_L(ISpin), C+1, C+A, 1, C )
Call GA_Scale_Patch(g_L(ISpin), C+1, C+A, 1, C, -1.0d0 )
C
C L(ay) = L2(ay) + L3(ay) + L4(ay)
C Total size ap pn pq
C
Call GA_Add_Patch(
$ 1.0d0, g_L2(ISpin), 1, NVir(ISpin), 1, NFrzO,
$ 1.0d0, g_L3(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), 1, NFrzO,
$ g_L(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), 1, NFrzO)
Call GA_Add_Patch( 1.0d0, g_L4(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), 1, NFrzO,
$ 1.0d0, g_L(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), 1, NFrzO,
$ g_L(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), 1, NFrzO)
Write (6, *) 'Add patch in L(ay)'
Call GA_Add_Patch( 1.0d0, g_L2(ISpin), 1, V, 1, C,
$ 1.0d0, g_L3(ISpin), C+A+1, C+A+V, 1, C,
$ g_L(ISpin), C+A+1, C+A+V, 1, C)
Call GA_Add_Patch( 1.0d0, g_L4(ISpin), C+A+1, C+A+V, 1, C,
$ 1.0d0, g_L(ISpin), C+A+1, C+A+V, 1, C,
$ g_L(ISpin), C+A+1, C+A+V, 1, C)
EndIf ! NFrzO .gt. 0
C
C L(ai) = L1(ai) + L2(ai) + L3(ai) + L4(ai)
C Total size pi ap pn pq
C
Call GA_Add_Patch(
$ 1.0d0, g_L1(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), 1, NAct(ISpin),
$ 1.0d0, g_L2(ISpin), 1, NVir(ISpin), NFrzO+1,
$ NFrzO+NAct(ISpin),
$ g_L(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), NFrzO+1, NFrzO+NAct(ISpin))
Call GA_Add_Patch(
$ 1.0d0, g_L3(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), NFrzO+1, NFrzO+NAct(ISpin),
$ 1.0d0, g_L(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), NFrzO+1, NFrzO+NAct(ISpin),
$ g_L(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), NFrzO+1, NFrzO+NAct(ISpin))
Call GA_Add_Patch(
$ 1.0d0, g_L4(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), NFrzO+1, NFrzO+NAct(ISpin),
$ 1.0d0, g_L(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), NFrzO+1, NFrzO+NAct(ISpin),
$ g_L(ISpin), NFrzO+NAct(ISpin)+1,
$ NFrzO+NAct(ISpin)+NVir(ISpin), NFrzO+1, NFrzO+NAct(ISpin))
Write (6, *) 'Add patch in L(ai)'
Call GA_Add_Patch( 1.0d0, g_L1(ISpin), C+A+1, C+A+V, 1, A,
$ 1.0d0, g_L2(ISpin), 1, V, C+1, C+A,
$ g_L(ISpin), C+A+1, C+A+V, C+1, C+A)
Call GA_Add_Patch( 1.0d0, g_L3(ISpin), C+A+1, C+A+V, C+1, C+A,
$ 1.0d0, g_L(ISpin), C+A+1, C+A+V, C+1, C+A,
$ g_L(ISpin), C+A+1, C+A+V, C+1, C+A)
Call GA_Add_Patch( 1.0d0, g_L4(ISpin), C+A+1, C+A+V, C+1, C+A,
$ 1.0d0, g_L(ISpin), C+A+1, C+A+V, C+1, C+A,
$ g_L(ISpin), C+A+1, C+A+V, C+1, C+A)
C
C Print results if requested
C
If ( Util_Print('partial l', Print_Debug) )
$ Call GA_Print( g_L(ISpin) )
C
EndDo ! ISpin
C
C Print results if requested
C ***********************
C * L contributes to P2 *
C ***********************
C This is separate from the CPHF equations
C
If ( Util_Print('partial l', Print_Debug) ) then
C P2(iy) = 1/2 L1(yi) (e(i) - e(y))^{-1}
C
C
If ( NFrzO .gt. 0) then
Do ISpin = 1, TopSpin
Call GA_Print( g_L(ISpin) )
C
C = NFrzO
A = NAct(ISpin)
V = NVir(ISpin)
C
Call GA_Copy_Patch('T', g_L1(ISpin), 1, C, C+1, C+A,
$ g_P2(ISpin), C+1, C+A, 1, C)
Call GA_Scale_Patch(g_P2(ISpin), C+1, C+A, 1, C, 0.5d0 )
C
If ( Util_Print('partial p2', Print_Debug) )
$ Call GA_Print( g_P2(ISpin) )
C
C Find out what portion of matrix we hold
C
Call GA_Distribution(g_P2(ISpin), Me, PLo, PHi, QLo, QHi)
C
C If we own part of the AC patch, we have work to do
C
PLo = Max(PLo, C+1)
PHi = Min(PHi, C+A)
QLo = Max(QLo, 1)
QHi = Min(QHi, C)
If ( (PHi-PLo+1) * (QHi-QLo+1) .gt. 0 ) then
Call GA_Access(g_P2(ISpin), PLo, PHi, QLo, QHi,
$ Index, LD)
Do Y = QLo, QHi
Do I = PLo, PHi
Scale = 1/( Eig(I, ISpin) - Eig(Y, ISpin) )
Write ( 6, *) 'I, Y, Scale = ', I, Y, Scale
Dbl_MB( Index ) = Dbl_MB(Index) * Scale
Index = Index + 1
EndDo
Index = Index + LD - 1
EndDo
Call GA_Release_Update(g_P2(ISpin), PLo, PHi, QLo, QHi)
EndIf
Call GA_Sync
C
C Print results if requested
C
If ( Util_Print('partial p2', Print_Debug) )
$ Call GA_Print( g_P2(ISpin) )
C
EndDo
EndIf
C
C ***********************
C * L contributes to W2 *
C ***********************
C
Do ISpin = 1, TopSpin
Call GA_Zero( g_W2(ISpin) )
C
C = NFrzO
A = NAct(ISpin)
V = NVir(ISpin)
C
C NOTE: All L contributions to W2 come with a 1/2 factor attached.
C Also, all L contributions to the (C+A)x(C+A) region carry a minus
C sign. These will be handled at the end, after the matrix is
C assembled.
C
C W2(mi) <-- L1(mi)
C
Write (6, *) 'Add patch in W2(mi)'
Call GA_Add_Patch(1.0d0, g_L1(ISpin), 1, C+A, 1, A,
$ 1.0d0, g_W2(ISpin), 1, C+A, C+1, C+A,
$ g_W2(ISpin), 1, C+A, C+1, C+A)
C
C W2(ap) <-- L2(ap)
C
Write (6, *) 'Add patch in W2(ap)'
Call GA_Add_Patch(1.0d0, g_L2(ISpin), 1, V, 1, C+A+V,
$ 1.0d0, g_W2(ISpin), C+A+1, C+A+V, 1, C+A+V,
$ g_W2(ISpin), C+A+1, C+A+V, 1, C+A+V)
C
C W2(mi) <-- L3(mi) (Note: L3(ij) = L3(ji))
C
Write (6, *) 'Add patch in W2(mi)'
Call GA_Add_Patch(1.0d0, g_L3(ISpin), 1, C+A, 1, A,
$ 1.0d0, g_W2(ISpin), 1, C+A, C+1, C+A,
$ g_W2(ISpin), 1, C+A, C+1, C+A)
C
C W2(yz) <-- L3(yz)
C
Write (6, *) 'Add patch in W2(yz)'
If ( C .gt. 0)
$ Call GA_Add_Patch(1.0d0, g_L3(ISpin), 1, C, 1, C,
$ 1.0d0, g_W2(ISpin), 1, C, 1, C,
$ g_W2(ISpin), 1, C, 1, C)
C
C W2(om) <-- L4(om)
C
Write (6, *) 'Add patch in W2(om)'
Call GA_Add_Patch(1.0d0, g_L4(ISpin), 1, C+A, 1, C+A,
$ 1.0d0, g_W2(ISpin), 1, C+A, 1, C+A,
$ g_W2(ISpin), 1, C+A, 1, C+A)
C
C Handle minus sign and 1/2 factor
C
Call GA_Scale_Patch(g_W2(ISpin), 1, C+A, 1, C+A, -1.0d0)
Call GA_Scale(g_W2(ISpin), -0.5d0)
C
C Now handle the two blocks which we know from symmetry arguments
C but are not currently correct in W2.
C
C W2(ma) = W2(am)
C
Call GA_Copy_Patch('T',g_W2(ISpin), C+A+1, C+A+V, 1, C+A,
$ g_W2(ISpin), 1, C+A, C+A+1, C+A+V)
C
C W2(iy) = W2(yi)
C
If ( C .gt. 0 ) Call GA_Copy_Patch('T',
$ g_W2(ISpin), 1, C, C+1, C+A,
$ g_W2(ISpin), C+1, C+A, 1, C)
C
C
C Print results if requested
C
If ( Util_Print('partial w2', Print_Debug) )
$ Call GA_Print( g_W2(ISpin) )
C
EndDo ! ISpin
C
C Clean up memory
C
Do ISpin = 1, TopSpin

12
src/rimp2_grad/rimp2g.F
View file

@ -1,5 +1,5 @@
Logical Function RIMP2G( RTDB )
C$Id: rimp2g.F,v 1.7 1998-07-14 16:18:37 bernhold Exp $
C$Id: rimp2g.F,v 1.8 1998-07-14 21:41:49 bernhold Exp $
Implicit NONE
Integer RTDB
C
@ -374,7 +374,7 @@ C Clean up some memory
C
Status = .TRUE.
Status = Status .AND. MA_Pop_Stack( H_Contrib)
Status = Status .AND. MA_Free_Heap(H_Eval)
c$$$ Status = Status .AND. MA_Free_Heap(H_Eval)
If ( .NOT. Status) Call ErrQuit(
$ 'RIMP2G: Unable to destroy local arrays', 0)
C
@ -384,7 +384,8 @@ C ************
C
If ( DoGrad ) then
Call RIMP2_Driver_G(MaxSpin, TopSpin, NFrzOcc, NAct, NVir,
$ FitBas, FN_Int, FN_Gam, BraKetSame)
$ FitBas, FN_Int, FN_Gam, BraKetSame, g_P2,
$ Dbl_MB(I_Eval), MxNCorBF)
EndIf
C
C ***********
@ -393,6 +394,11 @@ C ***********
C
If ( DRA_Terminate() .ne. 0) Call ErrQuit(
$ 'RIMP2G: DRA_Terminate failed', 0)
C
Status = .TRUE.
Status = Status .AND. MA_Free_Heap(H_Eval)
If ( .NOT. Status) Call ErrQuit(
$ 'RIMP2G: Unable to destroy local arrays', 0)
C
Status = .TRUE.
Status = Status .AND. Geom_Destroy( Geom)