forked from hphratchian/spinContaminationChecker
-
Notifications
You must be signed in to change notification settings - Fork 0
/
spinContaminationChecker1.f03
508 lines (495 loc) · 18.4 KB
/
spinContaminationChecker1.f03
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
INCLUDE 'hphSpinFun_mod.f03'
Program spinContaminationChecker
!
! This program evaluates the full set of S^2 eigenvalues and eigenvectors in
! the basis of substituted (Sz conserving) determinants. The program reads a
! matrix element file generated from a Gaussian unrestricted HF or DFT
! calculation.
!
! This code is based on a previous program written by Hratchian and relies
! on S2 matrix element code originally written by Lee Thompson.
!
!
! -H. P. Hratchian, 2021.
!
!
! USE Connections
!
use hphSpinFun_mod
USE OMP_LIB
!
! Variable Declarations
!
implicit none
integer(kind=int64)::nCommands,nSwitches,iPrint,nOMP,i,j,k,k1,k2, &
nAtoms,nAt3,nBasis,nDetAlpha,nDetBeta,nDetTotal,NDetTT, &
nBit_Ints,iAlpha,iBeta,jAlpha,jBeta,nFrozenCore,nFrozenVirtual
integer(kind=int64),dimension(:),allocatable::atomicNumbers, &
tmpVectorInt1,tmpVectorInt2,tmpVectorInt3,tmpVectorInt4
integer(kind=int64),dimension(:,:),allocatable::permuteAlpha, &
permuteBeta
integer(kind=int64),dimension(:,:,:),allocatable::stringLeftAlpha, &
stringLeftBeta,stringRightAlpha,stringRightBeta
real(kind=real64)::t1,t2,t1A,t2A,Vnn,Escf
real(kind=real64),dimension(3)::tmp3Vec
real(kind=real64),dimension(:),allocatable::cartesians, &
S2_MatSymm,tmpEVals,tmpVector
real(kind=real64),dimension(:,:),allocatable::distanceMatrix, &
S2_Mat,tmp2NBasisSq,tmpEVecs,tmpMatrix
character(len=512)::matrixFilename,tmpString
logical::DEBUG=.false.,gaussianCall,correspondingOrbitals, &
ABOverlapTest
type(mqc_gaussian_unformatted_matrix_file)::GMatrixFile
type(MQC_Variable)::tmpMQCvar
type(MQC_Scalar)::tmpMQCvar1,tmpMQCvar2,tmpMQCvar3,tmpMQCvar4, &
tmpMQCvar5
type(MQC_Variable)::nEalpha,nEbeta,nEtot,KEnergy,VEnergy,OneElEnergy, &
TwoElEnergy,scfEnergy,SSqTmp,SzTotal,MultiplicityTemp
type(MQC_Variable)::SMatrixAO,SMatrixMOAB,TMatrixAO,VMatrixAO, &
HCoreMatrixAO,FMatrixAlpha,FMatrixBeta,PMatrixAlpha, &
PMatrixBeta,PMatrixTotal,ERIs,JMatrix,KMatrixAlpha, &
CAlpha,CBeta
type(MQC_R4Tensor)::tmpR4
Type(MQC_Determinant)::Determinants
type(MQC_Variable)::mqcMat_pA,mqcMat_pB,mqcMat_tmp2NBS,mqcMat_S2, &
mqcMat_S2_Symm,mqcMat_tmpEVals,mqcMat_tmpEVecs, &
tempMQCvar1,tempMQCvar2
!
! Format Statements
!
1000 Format(1x,'Enter Test Program scfEnergyTerms.')
1010 Format(3x,'Matrix File: ',A)
1020 Format(3x,'Print Level: ',I2)
1030 Format(3x,'Number OMPs: ',I3,/)
1100 Format(1x,'nAtoms=',I4)
1200 Format(1x,'Atomic Coordinates (Angstrom)')
1210 Format(3x,I3,2x,A2,5x,F7.4,3x,F7.4,3x,F7.4)
1300 Format(1x,'Nuclear Repulsion Energy = ',F20.6)
2000 Format(1x,60('-'),/,1x,'left strings: (',I7,' | ',I7,')', &
2x,'==> (',B64,' | ',B64,' )',/, &
1x,'right strings: (',I7,' | ',I7,')',2x, &
'==> (',B64,' | ',B64,' )')
2010 Format(1x,'S^2 = ',F15.6)
2100 Format(1x,'Non-Zero Projection (',F10.4,') onto Reference by EVec ', &
I10,'. EVal = ',F10.4)
5000 Format(1x,A,' CPU Time: ',F12.2,' s.')
8999 Format(/,1x,'END OF TEST PROGRAM scfEnergyTerms.')
!
!
! Start the job timer and let the user know we are starting the program.
!
call cpu_time(t1)
write(IOut,1000)
call flush(iOut)
!
! Read the command line arguments.
!
nCommands = command_argument_count()
if(nCommands.lt.1) &
call mqc_error('No command line arguments provided. At least one command line argument is required giving the input Gaussian matrix element file name.')
!hph+
write(*,*)
write(*,*)
call commandLineArgs_switches(nSwitches,gaussianCall, &
correspondingOrbitals,ABOverlapTest,nFrozenCore,nFrozenVirtual, &
permuteAlpha,permuteBeta)
write(*,*)' Hrant - nSwitches = ',nSwitches
write(*,*)' gaussianCall = ',gaussianCall
write(*,*)' correspondingOrbitals = ',correspondingOrbitals
write(*,*)' nFrozenCore = ',nFrozenCore
write(*,*)' nFrozenVirtual = ',nFrozenVirtual
!CJD+
! mqcMat_pA = permuteAlpha
! call mqcMat_pA%print(header='CJD - Alpha Permutations')
! call mqc_print(iOut,mqcMat_pA,header='Alpha Permutations')
! mqcMat_pB = permuteBeta
! call mqcMat_pB%print(header='CJD - Beta Permutations')
! call mqc_print(iOut,mqcMat_pB,header='Beta Permutations')
!CJD-
!hph-
!hph+
! call get_command_argument(1,tmpString)
! if(TRIM(tmpString).eq.'-g') then
! gaussianCall = .true.
! else
! gaussianCall = .false.
! endIf
!hph-
if(gaussianCall) then
call commandLineArgs_gaussian(iOut,matrixFilename,iPrint,nOMP)
else
call commandLineArgs_direct(nSwitches,matrixFilename,iPrint,nOMP)
endIf
write(iOut,*)' matrixFilename: ',TRIM(matrixFilename)
write(iOut,1020) iPrint
if(iPrint.eq.-1) then
iPrint = 10
DEBUG = .true.
endIf
write(iOut,1030) nOMP
!hph+
! goto 999
!hph-
call omp_set_num_threads(nOMP)
!
! Open the Gaussian matrix file and load the number of atomic centers.
!
call GMatrixFile%load(matrixFilename)
write(IOut,1010) TRIM(matrixFilename)
nAtoms = GMatrixFile%getVal('nAtoms')
write(IOut,1100) nAtoms
!
! Figure out nAt3, then allocate memory for key arrays.
!
nAt3 = 3*nAtoms
Allocate(cartesians(NAt3),atomicNumbers(NAtoms))
!
! Set the intrinsic nBasis and allocate space for the temp intrinsic nBasis
! x nBasis matrix.
!
nBasis = Int(GMatrixFile%getVal('nbasis'))
!hph if(nBasis.gt.60) call mqc_error('NBasis > 60 NYI.')
Allocate(tmp2NBasisSq(2*nBasis,2*nBasis))
!
! Load up a few matrices from the matrix file.
!
call GMatrixFile%getArray('OVERLAP',mqcVarOut=SMatrixAO)
call GMatrixFile%getArray('KINETIC ENERGY',mqcVarOut=TMatrixAO)
call GMatrixFile%getArray('CORE HAMILTONIAN ALPHA',mqcVarOut=HCoreMatrixAO)
call GMatrixFile%getArray('ALPHA FOCK MATRIX',mqcVarOut=FMatrixAlpha)
if(GMatrixFile%isUnrestricted()) then
call GMatrixFile%getArray('BETA FOCK MATRIX',mqcVarOut=FMatrixBeta)
else
FMatrixBeta = FMatrixAlpha
endIf
call GMatrixFile%getArray('ALPHA SCF DENSITY MATRIX',mqcVarOut=PMatrixAlpha)
call GMatrixFile%getArray('ALPHA MO COEFFICIENTS',mqcVarOut=CAlpha)
if(GMatrixFile%isUnrestricted()) then
call GMatrixFile%getArray('BETA SCF DENSITY MATRIX',mqcVarOut=PMatrixBeta)
call GMatrixFile%getArray('BETA MO COEFFICIENTS',mqcVarOut=CBeta)
else
PMatrixBeta = PMatrixAlpha
CBeta = CAlpha
endIf
PMatrixTotal = PMatrixAlpha+PMatrixBeta
VMatrixAO = HCoreMatrixAO-TMatrixAO
!
! Make permutations for the alpha and beta MO coefficients.
!
write(*,*)
write(*,*)
call CAlpha%print(iOut,header='CAlpha BEFORE Permutations')
call CBeta%print(iOut,header='CBeta BEFORE Permutations')
if(Allocated(permuteAlpha)) then
do i = 1,SIZE(permuteAlpha,2)
call MQC_Variable_MatrixPermuteColumns(CAlpha, &
permuteAlpha(1,i),permuteAlpha(2,i))
endDo
endIf
if(Allocated(permuteBeta)) then
do i = 1,SIZE(permuteBeta,2)
call MQC_Variable_MatrixPermuteColumns(CBeta, &
permuteBeta(1,i),permuteBeta(2,i))
endDo
endIf
write(*,*)
write(*,*)
call CAlpha%print(iOut,header='CAlpha AFTER Permutations')
call CBeta%print(iOut,header='CBeta AFTER Permutations')
write(*,*)
write(*,*)
!hph goto 999
!
! Calculate the number of electrons using <PS>.
!
nEalpha = Contraction(PMatrixAlpha,SMatrixAO)
nEbeta = Contraction(PMatrixBeta,SMatrixAO)
nEtot = Contraction(PMatrixTotal,SMatrixAO)
nEalpha = GMatrixFile%getVal('nalpha')
nEbeta = GMatrixFile%getVal('nbeta')
call nEalpha%print(IOut,' <P(Alpha)S>=')
call nEbeta%print(IOut,' <P(Beta )S>=')
call nEtot%print(IOut,' <P(Total)S>=')
!
! Build the total energy.
!
scfEnergy = MQC(0.5)* &
(Contraction(PMatrixAlpha,(HCoreMatrixAO+FMatrixAlpha)) + &
Contraction(PMatrixBeta,(HCoreMatrixAO+FMatrixBeta)))
call scfEnergy%print(iOut,'electronic energy (no Vnn) = ')
!
! Calculate the 1-electron energy and component pieces of the 1-electron
! energy. Also, calculate the 2-electron energy.
!
KEnergy = Contraction(PMatrixTotal,TMatrixAO)
VEnergy = Contraction(PMatrixTotal,VMatrixAO)
OneElEnergy = Contraction(PMatrixTotal,HCoreMatrixAO)
TwoElEnergy = scfEnergy-OneElEnergy
call KEnergy%print(IOut,' <P.K> = ')
call VEnergy%print(IOut,' <P.V> =')
call OneElEnergy%print(IOut,' <P.H> = ')
call TwoElEnergy%print(IOut,' EE, <P.F>-<P.H> = ')
!
! Load the atommic numbers and Cartesian coordinates into our intrinsic
! arrays.
!
atomicNumbers = GMatrixFile%getAtomicNumbers()
cartesians = GMatrixFile%getAtomCarts()
cartesians = cartesians*angPBohr
!
! Print out the atomic numbers and Cartesian coordiantes for each atomic
! center.
!
write(IOut,1200)
do i = 1,NAtoms
j = 3*(i-1)
write(IOut,1210) i,mqc_element_symbol(atomicNumbers(i)), &
cartesians(j+1),cartesians(j+2),cartesians(j+3)
endDo
!
! Form the distance matrix between atomic centers.
!
Allocate(distanceMatrix(nAtoms,nAtoms))
do i = 1,nAtoms-1
distanceMatrix(i,i) = float(0)
k1 = 3*(i-1)+1
do j = i+1,NAtoms
k2 = 3*(j-1)+1
tmp3Vec = cartesians(k1:k1+2)-cartesians(k2:k2+2)
distanceMatrix(i,j) = sqrt(dot_product(tmp3Vec,tmp3Vec))
distanceMatrix(j,i) = distanceMatrix(i,j)
endDo
endDo
!
! Calculate the nuclear-nuclear repulsion energy.
!
distanceMatrix = distanceMatrix/angPBohr
Vnn = float(0)
do i = 1,NAtoms-1
do j = i+1,NAtoms
Vnn = Vnn + float(atomicNumbers(i)*atomicNumbers(j))/distanceMatrix(i,j)
endDo
endDo
write(iOut,1300) Vnn
!
! Put things together and report the SCF energy.
!
scfEnergy = scfEnergy + MQC(Vnn)
call scfEnergy%print(IOut,' Total SCF Energy = ')
call flush(iOut)
!
! Build the alpha/beta MO overlap matrix and print it out.
!
if(DEBUG) call SMatrixAO%print(header='AO Overlap Matrix')
SMatrixMOAB = MatMul(TRANSPOSE(CAlpha),MatMul(SMatrixAO,CAlpha))
if(DEBUG) call SMatrixMOAB%print(header='CAlpha.S.CAlpha')
SMatrixMOAB = MatMul(TRANSPOSE(CBeta),MatMul(SMatrixAO,CBeta))
if(DEBUG) call SMatrixMOAB%print(header='CBeta.S.CBeta')
SMatrixMOAB = MatMul(TRANSPOSE(CAlpha),MatMul(SMatrixAO,CBeta))
if(iPrint.ge.0.or.DEBUG) &
call SMatrixMOAB%print(header='Alpha-Beta MO Overlap Matrix')
if(ABOverlapTest) goto 999
!
! Try building things up for post-SCF like jobs.
!
write(iOut,*)
write(iOut,*)
write(iOut,*)' Hrant - About to call Gen_Det_Str.'
call flush(iOut)
tmpMQCvar1 = GMatrixFile%getVal('nbasis')
tmpMQCvar2 = INT(nEalpha)
tmpMQCvar3 = INT(nEbeta)
tmpMQCvar4 = nFrozenCore
tmpMQCvar5 = nFrozenVirtual
tmpMQCvar1 = tmpMQCvar1 - tmpMQCvar4 - tmpMQCvar5
tmpMQCvar2 = tmpMQCvar2 - tmpMQCvar4
tmpMQCvar3 = tmpMQCvar3 - tmpMQCvar4
call cpu_time(t1A)
Call Gen_Det_Str(IOut,2,tmpMQCvar1,tmpMQCvar2, &
tmpMQCvar3,Determinants,tmpMQCvar4)
call cpu_time(t2A)
write(iOut,5000) 'Gen_Det_Str',t2A-t1A
nDetAlpha = Bin_Coeff(GMatrixFile%getVal('nbasis')-nFrozenCore-nFrozenVirtual, &
INT(nEalpha)-nFrozenCore)
nDetBeta = Bin_Coeff(GMatrixFile%getVal('nbasis')-nFrozenCore-nFrozenVirtual, &
INT(nEbeta)-nFrozenCore)
nDetTotal = nDetAlpha*nDetBeta
write(iOut,*)
write(iOut,*)' nBasis = ',INT(GMatrixFile%getVal('nbasis'))
write(iOut,*)' nEalpha = ',INT(nEalpha)
write(iOut,*)' nEbeta = ',INT(nEbeta)
write(iOut,*)' nDetAlpha = ',nDetAlpha
write(iOut,*)' nDetBeta = ',nDetBeta
write(iOut,*)' nDetTotal = ',nDetTotal
write(iOut,*)
nBit_Ints = ((INT(GMatrixFile%getVal('nbasis'))-nFrozenVirtual)/(Bit_Size(0)-1))+1
write(iOut,*)' nBit_Ints = ',nBit_Ints
write(iOut,*)
call flush(iOut)
!hph+
! goto 999
!hph-
Allocate(S2_Mat(nDetTotal,nDetTotal))
tmp2NBasisSq(1:NBasis,1:NBasis) = float(0)
tmp2NBasisSq(NBasis+1:2*NBasis,NBasis+1:2*NBasis) = float(0)
call MQC_Variable_mqc2intrinsicReal2Array(tmpMatrix,SMatrixMOAB)
tmp2NBasisSq(1:NBasis,NBasis+1:2*NBasis) = tmpMatrix
tmp2NBasisSq(NBasis+1:2*NBasis,1:NBasis) = TRANSPOSE(tmpMatrix)
do i = 1,2*nBasis
tmp2NBasisSq(i,i) = float(1)
endDo
!CJD+
mqcMat_tmp2NBS = tmp2NBasisSq
if(iPrint.ge.0.or.DEBUG) &
call mqcMat_tmp2NBS%print(header='CJD - Full MO-MO Overlap')
! call mqc_print(iOut,mqcMat_tmp2NBS,header='Full MO-MO Overlap')
!CJD-
call flush(iOut)
!
! Pre-process the string list combinations for the S2 matrix element
! formation loops below.
!
call cpu_time(t1A)
Allocate(stringLeftAlpha(nDetTotal,nDetTotal,nBit_Ints), &
stringLeftBeta(nDetTotal,nDetTotal,nBit_Ints), &
stringRightAlpha(nDetTotal,nDetTotal,nBit_Ints), &
stringRightBeta(nDetTotal,nDetTotal,nBit_Ints), &
tmpVectorInt1(nBit_Ints),tmpVectorInt2(nBit_Ints), &
tmpVectorInt3(nBit_Ints),tmpVectorInt4(nBit_Ints))
call flush(iOut)
do iBeta = 1,NDetBeta
tmpVectorInt1 = Determinants%Strings%Beta%vat([iBeta], &
[1,nBit_Ints])
do iAlpha = 1,nDetAlpha
tmpVectorInt2 = Determinants%Strings%Alpha%vat([iAlpha], &
[1,nBit_Ints])
i = (iBeta-1)*nDetAlpha+iAlpha
stringLeftBeta(i,1,:) = tmpVectorInt1
stringLeftAlpha(i,1,:) = tmpVectorInt2
endDo
endDo
call flush(iOut)
do jBeta = 1,NDetBeta
tmpVectorInt3 = Determinants%Strings%Beta%vat([jBeta], &
[1,nBit_Ints])
do jAlpha = 1,nDetAlpha
j = (jBeta-1)*nDetAlpha+jAlpha
stringRightBeta(1,j,:) = tmpVectorInt3
tmpVectorInt4 = Determinants%Strings%Alpha%vat([jAlpha], &
[1,nBit_Ints])
stringRightAlpha(1,j,:) = tmpVectorInt4
endDo
endDo
call cpu_time(t2A)
write(iOut,5000) 'S2 Matrix String Pre-Processing',t2A-t1A
call flush(iOut)
!
! Evaluate the S2 matrix elements.
!
call cpu_time(t1A)
!$OMP PARALLEL DO DEFAULT(NONE), &
!$OMP SHARED(S2_Mat,stringLeftAlpha,stringLeftBeta,stringRightAlpha,stringRightBeta, &
!$OMP tmp2NBasisSq,NDetAlpha,NDetBeta,nBasis,iOut,DEBUG), &
!$OMP PRIVATE(i,j)
do i = 1,NDetAlpha*NDetBeta
! do j = 1,NDetAlpha*NDetBeta
do j = 1,i
S2_Mat(i,j) = S2_Mat_Elem1(IOut,2,nBasis, &
stringLeftAlpha(i,1,:),stringLeftBeta(i,1,:), &
stringRightAlpha(1,j,:),stringRightBeta(1,j,:),tmp2NBasisSq,1_int64)
! S2_Mat(i,j) = S2_Mat_Elem(IOut,2,nBasis, &
! stringLeftAlpha(i,1,:),stringLeftBeta(i,1,:), &
! stringRightAlpha(1,j,:),stringRightBeta(1,j,:),tmp2NBasisSq)
if(ABS(S2_Mat(i,j)).lt.(float(1)/float(10000))) S2_Mat(i,j) = float(0)
S2_Mat(j,i) = S2_Mat(i,j)
! if(DEBUG) write(iOut,2010) S2_Mat(i,j)
endDo
endDo
!$OMP END PARALLEL DO
call cpu_time(t2A)
write(iOut,5000) 'S2 Matrix Formation',t2A-t1A
DeAllocate(stringLeftAlpha,stringLeftBeta,stringRightAlpha, &
stringRightBeta)
!CJD+
mqcMat_S2 = S2_Mat
if(iPrint.ge.2.or.DEBUG) &
call mqcMat_S2%print(header='CJD - S2 Matrix (FULL)')
! call mqc_print(iOut,mqcMat_S2,header='S2 Matrix (FULL)')
!CJD-
call flush(iOut)
!
! Diagonalize S2_Mat and report the eigenvalues.
!
nDetTT = (nDetTotal*(nDetTotal+1))/2
Allocate(S2_MatSymm(nDetTT),tmpEVals(nDetTotal), &
tmpEVecs(nDetTotal,nDetTotal),tmpVector(3*nDetTotal))
k = 0
do i = 1,nDetTotal
do j = i,nDetTotal
k = k+1
S2_MatSymm(k) = S2_Mat(i,j)
endDo
endDo
!CJD+
mqcMat_S2_Symm = mqc_matrixSymm2Full(S2_MatSymm,'L')
if(DEBUG) &
call mqcMat_S2_Symm%print(header='CJD - S2 Matrix from Symm')
! call mqc_print(iOut,mqcMat_S2_Symm,header='S2 Matrix from Symm')
!CJD-
write(iOut,*)' k = ',k
write(iOut,*)' nDetTT = ',nDetTT
write(iOut,*)
write(iOut,*)
!hph goto 999
write(iOut,*)' Calling DSPEV.'
call cpu_time(t1A)
Call DSPEV('V','L',nDetTotal,S2_MatSymm,tmpEVals,tmpEVecs, &
NDetTotal,tmpVector,i)
call cpu_time(t2A)
write(iOut,5000) 'S2 Diagonalization',t2A-t1A
write(iOut,*)' After DSPEV, i = ',i
call flush(iOut)
do i = 1,nDetTotal
if(ABS(tmpEVals(i)).lt.(float(1)/float(10000))) tmpEVals(i) = float(0)
endDo
!CJD+
mqcMat_tmpEVals = tmpEVals
call mqcMat_tmpEVals%print(header='CJD - S2 EVals')
! call mqc_print(iOut,mqcMat_tmpEVals,header='S2 EVals')
mqcMat_tmpEVecs = tmpEVecs
if(iPrint.ge.2.or.DEBUG) &
call mqcMat_tmpEVecs%print(header='CJD - S2 EVecs')
! call mqc_print(iOut,mqcMat_tmpEVecs,header='S2 EVecs')
tempMQCvar1 = MatMul(Transpose(tmpEVecs),MatMul(S2_Mat,tmpEVecs))
if(DEBUG) &
call tempMQCvar1%print(header='CJD - EVecs.S2.EVecs')
! call mqc_print(iOut,tempMQCvar1,header='EVecs.S2.EVecs')
tempMQCvar2 = MatMul(Transpose(tmpEVecs),tmpEVecs)
if(DEBUG) &
call tempMQCvar2%print(header='CJD - EVecs.EVecs')
! call mqc_print(iOut,tempMQCvar2,header='EVecs.EVecs')
!CJD-
write(iOut,*)
write(iOut,*)' nDetTotal = ',nDetTotal
call flush(iOut)
do i = 1,NDetTotal
if(ABS(tmpEVecs(nDetTotal,i)).gt.(float(5)/float(100))) then
write(iOut,2100) tmpEVecs(nDetTotal,i),i,tmpEvals(i)
endIf
endDo
!
! Try to figure out the range of Sz values found in the S^2 eigenvalue
! range.
!
SzTotal = MQC(0.5)*(nEAlpha-nEBeta)
MultiplicityTemp = nEAlpha-nEBeta+MQC(1)
call SzTotal%print(header='Sz=')
call MultiplicityTemp%print(header='Multiplicity=')
call sumOverSpinStates(iOut,INT(MultiplicityTemp),tmpEVals, &
tmpEVecs(nDetTotal,:))
!
999 Continue
write(iOut,8999)
call cpu_time(t2)
write(iOut,5000) 'TOTAL JOB',t2-t1
end program spinContaminationChecker