-
Notifications
You must be signed in to change notification settings - Fork 0
/
RSPT_L.f90
2453 lines (2177 loc) · 114 KB
/
RSPT_L.f90
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
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!****************************************************************************
!
! PROGRAM: RSPT_L
!
! PURPOSE: Entry point for the console application.
!
! RSPT_L.f90
!
! FUNCTIONS:
! RSPT_L - Entry point of console application.
!
!****************************************************************************
PROGRAM RSPT_L
! USE MKL_SCALAPACK
IMPLICIT REAL*16 (A-H,O-Z), INTEGER*4 (I-N)
INTEGER*4 LOUT,NLOG,NINP,NOUT, NAUX,NCSR,NTMP
COMMON /IOLUN/ LOUT,NLOG,NINP,NOUT, NAUX,NSCR,NTMP
! INTEGER*4 :: NQ
REAL*8 :: RE, IM !
COMPLEX*16 ELEMAT, LZELE, CZERO, CONE, TERM, SCPRO, SUMVIR, SUMORD !
REAL *16 MATELE, MATELEHARM, MATELEPERT !
DIMENSION NHORD(0:4) ! Number of Terms in Hamiltonian (i)
DIMENSION COEFLADD(0:63,0:51) !
COMMON /LADCOE/ COEFLADD !
ALLOCATABLE :: HCOEF(:), VCOEF(:), HOCOEF(:) !
COMPLEX*16, ALLOCATABLE :: LCOEF(:), LZMAT(:,:) !
COMPLEX*16, ALLOCATABLE :: WORK(:), TRIDIAG(:,:), TAU(:), EIGVECLZ(:,:), AUX(:,:) !
INTEGER *8, ALLOCATABLE :: IWORK(:) !
REAL *8, ALLOCATABLE :: DIAG(:), RWORK(:), EIGVAL(:) , OFFDIAG(:), EIGVALLZ(:) !
INTEGER *1, ALLOCATABLE :: LADHAM(:,:,:), LADHO(:,:,:), LADVP(:,:,:) !
INTEGER *1, ALLOCATABLE :: LADLZ (:,:,:) !
INTEGER *4, ALLOCATABLE :: ISTATE(:) !
INTEGER *4, ALLOCATABLE :: IQUANT(:,:), JQUANT(:,:), IEXQUANT(:,:), JEXQUANT(:,:),SELEQUANT(:,:) !
INTEGER *4, ALLOCATABLE :: IOPBRA(:), IOPKET(:) !
INTEGER *4, ALLOCATABLE :: IVIRTU(:,:), NQUPOL(:) !
INTEGER *4, ALLOCATABLE :: VMAXS(:), AUX1(:,:), AUX2(:,:) !
COMPLEX*16, ALLOCATABLE :: HARM(:,:), VPERT(:,:), HARMTRAN(:,:), VPERTRAN(:,:) !
COMPLEX*16, ALLOCATABLE :: PROJECT(:,:), PROJTRAN(:,:), COPY1(:,:), COPY2(:,:), HAMATRE(:,:) !
COMPLEX*16, ALLOCATABLE :: AUX3(:,:), AUX4(:,:) !
REAL *16, ALLOCATABLE :: W(:), WEX(:),ENERHO(:) !
ALLOCATABLE KPOLY(:), KEXPOLY(:) !
LOGICAL, ALLOCATABLE :: LOGICSTATE(:)
INTEGER *4, ALLOCATABLE :: IASSIG(:), IDEG(:)
REAL *16, ALLOCATABLE :: ENERHARM(:), ENERCI(:), ENERPT(:) !
REAL *8, ALLOCATABLE :: ENERDEG(:) !
REAL *16, ALLOCATABLE :: EPTORD(:)
COMPLEX*16, ALLOCATABLE :: PSIFUN(:,:), VPERTDEG(:,:)
INTEGER*8 LSPLIT ! Overall time for all states
INTEGER*8 LSPLIT0 ! Overall time for a single state
INTEGER*8 LSPLIT1 ! VCI calculation for a single state
INTEGER*8 LSPLIT2 ! RSPT calculation for a single state
CHARACTER*16 TIMTXT
CHARACTER*80 :: FILE_H(0:4),FILE_LZ
CHARACTER*80 :: FILE_HTOT, FILE_POLY, FILE_POLY2, FILE_HARM, FILE_PERT,FILE_TEMP
CHARACTER*80 :: QUANID
CHARACTER*64 PATH_CUR,PATH_PRO,PATH_OUT,PATH_BIN
COMMON /PATHS/ PATH_CUR,PATH_PRO,PATH_OUT,PATH_BIN
CHARACTER*112 LINEWIDE
CHARACTER* 64 FILE_CFG,FILE_MOL,FILE_INP,FILE_LST,FILE_STATE
CHARACTER* 64 LINE,FILE,PATH
CHARACTER* 16 LPES(4)
CHARACTER* 8 KEYWOR
CHARACTER* 16 KEYWOR16
CHARACTER* 64 COMPUTER
CHARACTER*300 LINEMAT
CHARACTER* 6 SDTQ(0:4)
CHARACTER*24 DATVER, MOLECULE
CHARACTER* 2 VERTICL
LOGICAL :: FLAG,STATUS,SELECT
LOGICAL :: CONVER
DATA NLOG,NINP,NOUT,NAUX,NSCR,NTMP /1,2,3,4,7,8/
DATA LPES /' ', '(Quartic PES)', ' ', '(Sextic PES)'/
DATA SDTQ /'Polyad','Single','Double','Triple','Quadru'/
DATA VERTICL /'|'/
DATA LOUT /2/
DATA NOPER /2/
DATA TOLCUT /1.0D-6/
DATA TOLOUT /1.0D-8/ ! Convergence criterion
DATA LIMIT /300/
DATA ZERO,HALF,ONE,TWO,THREE,FOUR /0.0D0,0.5D0,1.0D0,2.0D0,3.0D0,4.0D0/
DATA CZERO, CONE /(0.0D0,0.0D0), (1.0D0,0.0D0)/
DATVER = 'Jul 17, 2021, Sat 9:00'
NUMVER = 0100
!MOLECULE = 'C2H2'!'2DIsoHO' !'C2H2' !'CO2' ! Type of Molecule
!NQ = 7 ! 4 ! Number of Normal Mode , 4 For CO2, 7 For C2H2
!NSAYVETZ = 0 ! Sayvetz Condition, L = 0, 1, 2,...
!ALLOCATE ( W (NQ ))
!W(1) = 100.0D0
!W(2) = 100.0D0
!W(3) = 110.0D0
!W(4) = 110.0D0
!W(1) = 1.0D0
!W(2) = 1.0D0
!W(1) = 0.1353779746D+04 ! Normal Frequencies of CO2
!W(2) = 0.2396532220D+04 !
!W(3) = 0.6728925740D+03 !
!W(4) = 0.6728925740D+03 !
!W(1) = 0.35068570D+04 ! Normal Frequencies of C2H2
!W(2) = 0.20111680D+04 !
!W(3) = 0.34143540D+04 !
!W(4) = 0.62168700D+03 !
!W(5) = 0.62168700D+03 !
!W(6) = 0.74861800D+03 !
!W(7) = 0.74861800D+03 !
PATH_CUR = 'C:/PROG/RSPT_L/'
FILE_CFG = 'RSPT_L.ini'
!-----------------------------------------------------------------------------------------------
! Read RSPT configuration file RSPT.dat:
! 1) Look for it in Current directory, otherwise in
! 2) C:/PROG/
!
! RSPT working directory = C:/PROG/RSPT/
! ANCO import data file = RSPT_H2O_MP2_cc-pVTZ.mol
! RSPT input parameters file = RSPT_H2O_MP2_cc-pVTZ.inp
!-----------------------------------------------------------------------------------------------
OPEN (UNIT = NAUX, ACTION = 'READ', FORM = 'FORMATTED', FILE = FILE_CFG, ERR = 1000)
GO TO 1010
1000 CONTINUE
FILE_CFG = TRIM(PATH_CUR) // TRIM(FILE_CFG)
WRITE (*,1020) FILE_CFG
OPEN (UNIT = NAUX, ACTION = 'READ', FORM = 'FORMATTED', FILE = FILE_CFG, ERR = 910)
1010 CONTINUE
1020 FORMAT (/'RSPTFM: Configuration file RSPT.dat not found in ', &
& 'current directory,'/8X,'attempting alternative directory: ',A)
DO WHILE (.NOT. EOF (NAUX))
READ (NAUX,"(A16,32X,A64)") KEYWOR16,LINE
IF (KEYWOR16(1:1) .EQ. '!' .OR. KEYWOR16 .EQ. ' ') CYCLE
IF (KEYWOR16(1:10) .EQ. '# COMPUTER') COMPUTER = LINE ! Computer specifications
IF (KEYWOR16(1:10) .EQ. '# PATH_CUR') PATH_CUR = LINE ! Read RSPT working directory
IF (KEYWOR16(1:10) .EQ. '# PATH_OUT') PATH_OUT = LINE ! Taylor series output .dat directory
IF (KEYWOR16(1:10) .EQ. '# PATH_BIN') PATH_BIN = LINE ! Working directory for binary data (operators), .bin
IF (KEYWOR16(1:10) .EQ. '# PARA_MOL') FILE_MOL = LINE ! ANCO import data path for *.mol
IF (KEYWOR16(1:10) .EQ. '# RSPT_INP') FILE_INP = LINE ! RSPT input parameters file, RSPT_mol_MP2_cc-pVTZ.inp
IF (KEYWOR16(1:10) .EQ. '# RSPT_LST') FILE_LST = LINE ! RSPT results output file, RSPT_CO2.txt
IF (KEYWOR16(1:10) .EQ. '# LIST_LVL') &
& READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) LOUT ! Output level (0 = min, 1 = regular, 2 = full)
IF (KEYWOR16(1:10) .EQ. '# ACCURACY') &
& READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) NACCURMP ! FM accuracy
ENDDO
REWIND (NAUX)
OPEN (UNIT = NOUT, ACTION = 'WRITE', FORM = 'FORMATTED', FILE = FILE_LST, ERR = 930)
WRITE (NOUT,1030) FILE_CFG
1030 FORMAT (/'[RSPT_L]'//'Program configuration file: ',A64/112('-'))
DO WHILE (.NOT. EOF (NAUX))
READ (NAUX,"(A)") LINEWIDE
WRITE (NOUT,"(A)") LINEWIDE
END DO
CLOSE (NAUX)
WRITE (NOUT,"(112('-'))")
!-----------------------------------------------------------------------------------------------
! Read parameters in .mol file
!-----------------------------------------------------------------------------------------------
FILE_MOL = TRIM(PATH_CUR) // TRIM(FILE_MOL)
OPEN (UNIT = NAUX, ACTION = 'READ', FORM = 'FORMATTED', FILE = FILE_MOL, ERR = 920)
DO WHILE (.NOT. EOF(NAUX))
READ (NAUX,1200) KEYWOR, LINE
IF (KEYWOR(1:1) .EQ. '!') CYCLE
IF (KEYWOR .EQ. 'MOLE_NAM') THEN
MOLECULE = TRIM(LINE)
ELSE IF (KEYWOR .EQ. 'SAYV_CON') THEN
READ (UNIT = LINE, FMT = '(I4)') NSAYVETZ
ELSE IF (KEYWOR .EQ. 'NORM_MOD') THEN
READ (UNIT = LINE, FMT = '(I4)') NQ
ALLOCATE(W(NQ))
DO I = 1, NQ
READ (NAUX, '(8X,F24.16)') W(I)
ENDDO
ELSE IF (KEYWOR .EQ. 'END_FILE') THEN
EXIT
ENDIF
ENDDO
CLOSE (NAUX)
!-----------------------------------------------------------------------------------------------
! Check for existance of the input file:
!-----------------------------------------------------------------------------------------------
FILE_INP = TRIM(PATH_CUR) // TRIM(FILE_INP)
INQUIRE (FILE = FILE_INP, EXIST = STATUS)
IF (STATUS) THEN
OPEN (UNIT = NINP, ACTION = 'READ', FORM = 'FORMATTED', FILE = FILE_INP, ERR = 940)
ELSE
WRITE (*,1100) FILE_INP
PAUSE 'RSPT: Error -- acknowledge before STOP: <ENTER>'
STOP
ENDIF
1100 FORMAT (/'Input data file: ',A/'-- not found, termination.')
!-----------------------------------------------------------------------------------------------
! Input control parameters:
! LUN = NINP can be available or not available
!-----------------------------------------------------------------------------------------------
ALLOCATE ( KPOLY(NQ) )
TOLCUT = 1.0D-06 ! Tolerance for CVPT operator terms, 1/Cm
TOLOUT = 1.0D-16 ! RSPT convergence criterion
MPES = 2 ! Quartic or Sextic force field, 2/4
MAX_PT = 210 ! Maximum order of RSPT
MAXVRT = 10 ! Maximum mode excitation of virtual states
MULMOD = 3 ! Maximum modes involved in virtual states
ENERMAX = 20000.0 ! Energy limit for a virtual state
EFREMAX = 1000.0 ! Energy limit for reference states
NOPTST = 4 ! Maximum total quanta of excitation
NSTART = 0 ! Maximum total quanta of excitation
MAXPOL = 0 ! Polyad structure is undefined
!---- IF (MAXQUA .EQ. 0) MAXQUA = MAXPOL
MAXQUA = 0 ! Total quanta of excitation in a polyad
KPOLY = 0 ! Polyad coefficients initialization (Array)
NQEX = 0 ! The Number of External(NQEX, v_i) Quantum Number
NQIN = 0 ! The Number of Internal(NQIN, l_i) Quantum Number
INQUIRE (UNIT = NINP, OPENED = STATUS)
IF (STATUS) THEN
DO WHILE (.NOT. EOF (NINP))
READ (NINP,1200) KEYWOR,LINE
IF (KEYWOR(1:1) .EQ. '!' .OR. KEYWOR .EQ. ' ') CYCLE
IF (KEYWOR(1:4) .EQ. 'PARA') EXIT ! READ: 'PARA'
END DO
DO WHILE (.NOT. EOF (NINP))
READ (NINP,1200) KEYWOR,LINE
IF (KEYWOR(1:1) .EQ. '!' .OR. KEYWOR .EQ. ' ') CYCLE
IF (KEYWOR .EQ. 'CVPT_TOL') THEN
READ (UNIT = LINE, FMT = "(F8.0)", IOSTAT = IOS) TOLCUT
ELSE IF (KEYWOR .EQ. 'CONV_TOL') THEN
READ (UNIT = LINE, FMT = "(F8.0)", IOSTAT = IOS) TOLOUT
ELSE IF (KEYWOR .EQ. 'FORC_ORD') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) MPES
ELSE IF (KEYWOR .EQ. 'STAT_OPT') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) NOPTST
ELSE IF (KEYWOR .EQ. 'STAT_BEG') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) NSTART
ELSE IF (KEYWOR .EQ. 'RSPT_MAX') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) MAX_PT
ELSE IF (KEYWOR .EQ. 'VIRT_MAX') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) MAXVRT
ELSE IF (KEYWOR .EQ. 'MULT_MOD') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) MULMOD
ELSE IF (KEYWOR .EQ. 'ENER_MAX') THEN
READ (UNIT = LINE, FMT = "(F8.0)", IOSTAT = IOS) ENERMAX
ELSE IF (KEYWOR .EQ. 'EREF_MAX') THEN
READ (UNIT = LINE, FMT = "(F8.0)", IOSTAT = IOS) EREFMAX
ELSE IF (KEYWOR .EQ. 'EXTR_QUA') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) NQEX
ALLOCATE( KEXPOLY( NQEX ) )
ELSE IF (KEYWOR .EQ. 'INTR_QUA') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) NQIN
ELSE IF (KEYWOR .EQ. 'POLY_MAX') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) MAXPOL
IF (MAXPOL .GT. 0) THEN
READ (NINP,"(15I4)") (KPOLY (I),I=1,NQ)
READ (NINP,"(15I4)") (KEXPOLY(I),I=1,NQEX)
ENDIF
ELSE IF (KEYWOR .EQ. 'POLY_QUA') THEN
READ (UNIT = LINE, FMT = "(I4)", IOSTAT = IOS) MAXQUA
ENDIF
IF (KEYWOR(1:4) .EQ. 'VIBR' .OR. &
& KEYWOR(1:4) .EQ. 'DATA') THEN
BACKSPACE (NINP)
EXIT
ENDIF
IF (KEYWOR(1:4) .EQ. 'TERM') THEN ! Carry on reading input
EXIT
ENDIF
IF (KEYWOR(1:4) .EQ. 'STOP') THEN ! Close input file
CLOSE (NINP)
EXIT
ENDIF
END DO
ELSE
WRITE (*,1210)
ENDIF
!---- LUN = NINP stopped at STOP/EXIT or before VIBR/DATA.
1200 FORMAT (A8,3X,A64) ! 1400
1210 FORMAT (/'RSPT: Input not found, using default parameters.')
MULMOD = MIN (NQ,MULMOD) !! Attention should be paid in MIN
!-----------------------------------------------------------------------------------------------
! Options for selected states calculation
!-----------------------------------------------------------------------------------------------
INQUIRE (UNIT = NINP, OPENED = STATUS)
IF (STATUS) THEN
READ (NINP, "(A4)", ERR = 270) CONTROL
IF (CONTROL .EQ. 'VIBR') THEN
IF (NOPTST .GE. 0) STOP 'NSTOPT = -1 expected for VIBR'
NOPTST = 0 ! Overrides standard user-defined setting
ALLOCATE (SELEQUANT(NQ,0:255), STAT = IERR)
LOCERR = 3; IF (IERR .NE. 0) GO TO 990
DO Q = 1, NQ
SELEQUANT(Q,0) = 0
END DO
DO I = 1, 255
READ (NINP,"(15I4)", ERR = 260) (SELEQUANT(Q,I),Q=1,NQ)
NSELSTATE = I
END DO
!---- Input ends by appearing of next Keyword
260 BACKSPACE (NINP)
WRITE (NOUT,2110) NSELSTATE
DO I = 0, NSELSTATE
WRITE (NOUT,2120) I,(SELEQUANT(Q,I),Q=1,NQ)
END DO
ELSE
WRITE (*,2130) NOPTST
BACKSPACE (NINP)
ENDIF
ENDIF
270 CONTINUE
2110 FORMAT (/'Total number of user-defined states =',I4)
2120 FORMAT ('State number (',I4,') -- Quantum numbers:'/(5I4,2X,5I4))
2130 FORMAT (/'Specific vibrational states not provided, ', & ! 2100
& 'using all states up to ',I2,' quanta.')
!----------------------------------------------------------------------------
WRITE (NOUT,1300) REAL(NUMVER)/100.0,DATVER,MOLECULE
CALL TIMDAT (NOUT)
WRITE (*,1300) REAL(NUMVER)/100.0,DATVER,MOLECULE
CALL TIMDAT (0)
! INTENS = IOPT(3) ! >>>
! ISCATT = IOPT(4) ! >>>
1300 FORMAT (/72('=')/ &
& 'Large-Order Vibrational Rayleigh-Schroedinger Perturbation ', &
& 'Theory'/ &
& 'for expanding energies and wave functions into Taylor ', &
& 'series: Linear molecule'/72('=')// &
& 'Version : ',F6.2/ &
& 'Date : ',A24/ &
& 'Molecule: ',A64/)
!---- Initialize filenames for different orders of the Hamiltonian
NHORD = 0 ! Vector assignment
DO N = 0, MPES
WRITE (FMT = "('HAMILT_',A,'_H',I1,'.txt')", UNIT = FILE) TRIM(MOLECULE), N
FILE_H(N) = TRIM(PATH_BIN)//TRIM(FILE)
END DO
WRITE (NOUT,1310) ! Operators of Creation/Annihilation info.
WRITE (NOUT,1320) ! Operators of Creation/Annihilation info.
WRITE (FMT = "('H_TOTAL.BIN')", UNIT = FILE)
FILE_HTOT = TRIM(PATH_BIN)// TRIM(FILE)
1310 FORMAT (/'Second Quantization Hamiltonian will be input from ', &
& 'Wolfram Mathematica output File' )
1320 FORMAT (/'Express Hamiltonian in Creation/Annihilation ', &
& 'operators:'// &
& '(a) Creation of Quantum: a(k) = (q(k)-i*p(k))/2^(1/2)'/ &
& '(b) Annihilation of Quantum: b(k) = (q(k)+i*p(k))/2^(1/2)')
DO I = 0, MPES
OPEN(UNIT=NINP, FILE=FILE_H(I), FORM = "FORMATTED", ACTION = "READ")
NTERH = 0
DO WHILE (.NOT. EOF(NINP))
READ(NINP, "(I4,2X,F15.8)", IOSTAT = IOS) N, COEF
NTERH = NTERH + 1
END DO
NHORD(I) = NTERH
ENDDO
NTHAM = NHORD(0) + NHORD(1) + NHORD(2) + NHORD(3) + NHORD(4)
NHCOR = 0
WRITE (NOUT,1330) (NHORD(I),I=0,2),NHCOR,(NHORD(I),I=3,4)
1330 FORMAT (/'The number of operator terms in the Hamiltonian:--'/ &
& 'Zero Order (Harmonic) = ',I8/ &
& 'First Order (Cubic) = ',I8/ &
& 'Second Order (Quartic) = ',I8/ &
& 'Second Order (Coriolis) = ',I8/ &
& 'Third Order (Quintic) = ',I8/ &
& 'Fourth Order (Sextic) = ',I8)
!!---------------------------------------------------------------------------!!
! !!
!!---- Assemble all Orders of the Hamiltonian in FILE_HTOT[1..NTHAM] !!
! CALL COPYPOL(NQ, NHORD(0), FILE_H(0), FILE_HTOT) !!
!!---- NOTE: Routine ADDPOLF also collects like terms !!
! CALL ADDPOLF(NQ, NTHAM, FILE_HTOT, TOLCUT, NTHAM, FILE_HTOT, & !!
! & ONE, NHORD(1), FILE_H(1), ONE) !!
! CALL ADDPOLF(NQ, NTHAM, FILE_HTOT, TOLCUT, NTHAM, FILE_HTOT, & !!
! & ONE, NHORD(2), FILE_H(2), ONE) !!
!!---- Append Coriolis Term ! It will be done later !!
!!---- Append Sextic Terms ! It will be done later !!
! WRITE (NOUT,1340) MPES,NTHAM !!
!1340 FORMAT (/80('-')// & !!
! & 'The Vibrational Watson Hamiltonian (J = 0) is composed:'// & !!
! & 'The order of potential function =',I4/ & !!
! & 'Total number of primitive terms =',I12) !!
! !!
! !!
!!---- Allocate Memory and Upload the Hamiltonian !!
! ALLOCATE (HCOEF(NTHAM),LADHAM(NQ,NOPER,NTHAM), STAT = IERR) !!
! LOCERR = 2; IF (IERR .NE. 0) GO TO 990 !!
! !!
! !!
! CALL UPLOAD (NQ,NTHAM,FILE_HTOT,HCOEF,LADHAM,TOLCUT,LIMIT,0) !!
!!---- Arrange operator terms in the descending order of coefficients !!
! CALL SORTPOL (NQ,NOPER,NTHAM,HCOEF,LADHAM,TOLCUT,LIMIT,IERR) !!
! IF (IERR .GT. 0) STOP '[SORTPOL] Returned error code.' !!
! CALL OPEROUT (NQ,NOPER,NTHAM,HCOEF,LADHAM,TOLCUT,LIMIT) !!
!----------------------------------------------------------------------------!!
!---- Zero-order Hamiltonian
ALLOCATE (HOCOEF(NHORD(0)), LADHO(NQ,NOPER,NHORD(0)), STAT = IERR)
LOCERR = 14; IF (IERR .NE. 0) GO TO 990
FILE_HARM = TRIM(PATH_BIN) // 'H_HARM.BIN'
OPEN (UNIT = NAUX, FILE = FILE_HARM, STATUS = 'REPLACE', &
& ACTION = 'WRITE', BUFFERED = 'YES', FORM = 'UNFORMATTED')
CALL COPYPOL(NQ, NHORD(0), FILE_H(0), FILE_HARM)
!---- UPLOAD () calls OPEROUT () if IOUT = 1:
CALL UPLOAD (NQ,NHORD(0),FILE_HARM,HOCOEF,LADHO,TOLCUT,LIMIT,1) ! UPLOAD ()
CALL SORTPOL (NQ,NOPER,NHORD(0),HOCOEF,LADHO,TOLCUT,LIMIT,IERR) ! SORTPOL ()
IF (IERR .GT. 0) STOP '[ARRANGE] Returned error code.'
! WRITE(NOUT, )
!---- Pertubated Hamiltonian
NPERT = NHORD(1) + NHORD(2) ! + NHORD(3) + NHORD(4)
ALLOCATE (VCOEF(NPERT), LADVP(NQ,NOPER,NPERT), STAT = IERR)
LOCERR = 15; IF (IERR .NE. 0) GO TO 990
FILE_PERT = TRIM(PATH_BIN) // 'H_PERT.BIN'
OPEN (UNIT = NAUX, FILE = FILE_PERT, STATUS = 'REPLACE', &
& ACTION = 'WRITE', BUFFERED = 'YES', FORM = 'UNFORMATTED')
!----------------------------------------------------------------------------!!
! CALL COPYPOL(NQ, NTHAN, FILE_HTOT, FILE_PERT) !!
! CALL ADDPOLF(NQ, NPERT, FILE_PERT, TOLCUT, & !!
! & NTHAM, FILE_HTOT, ONE, NHORD(0), FILE_H(0), -ONE) !!
!----------------------------------------------------------------------------!!
CALL COPYPOL(NQ, NHORD(1), FILE_H(1), FILE_PERT)
CALL ADDPOLF(NQ, NPERT, FILE_PERT, TOLCUT, &
& NHORD(1), FILE_PERT, ONE, NHORD(2), FILE_H(2), ONE)
!---- UPLOAD () calls OPEROUT () if IOUT = 1:
CALL UPLOAD (NQ,NPERT,FILE_PERT, VCOEF,LADVP,TOLCUT,LIMIT,1) ! UPLOAD ()
CALL SORTPOL (NQ,NOPER,NERPT,VCOEF,LADVP,TOLCUT,LIMIT,IERR) ! SORTPOL ()
IF (IERR .GT. 0) STOP '[ARRANGE] Returned error code.'
!-------------------------------------------------------------------------------
! Prepare Angular Momentum L(z)
! Read Lz file from Mathematica Output File
!-------------------------------------------------------------------------------
WRITE (FMT = "('Lz_',A,'.txt')", UNIT = FILE) TRIM(MOLECULE)
FILE_LZ = TRIM(PATH_BIN)// TRIM(FILE)
OPEN(UNIT=NINP, FILE=FILE_LZ, FORM = "FORMATTED", ACTION = "READ")
NTERL = 0
DO WHILE (.NOT. EOF(NINP))
READ(NINP, "(I4,2(1X,F32.16))", IOSTAT = IOS) N, RE, IM
NTERL = NTERL + 1
END DO
REWIND(NINP)
NTLZ = NTERL
ALLOCATE(LCOEF (NTLZ ),LADLZ(NQ,2,NTLZ ))
NTERL = 1
DO WHILE (.NOT. EOF(NINP))
READ(NINP, "(I4,2(1X,F32.16),20(1X,I1))", IOSTAT = IOS) N, RE, IM, &
& ((LADLZ(Q,J,NTERL),J=1,NOPER),Q=1,NQ)
LCOEF(NTERL) = CMPLX(RE,IM, KIND = 8)
NTERL = NTERL + 1
ENDDO
!-----------------------------------------------------------------------------------------------
! Create/read specifications of vibrational states of interest: ! 2000
! IQUANT(vector, j = 0..NSTATE)
! j = 0: ground state
! j > 0: j(polyad), polyad = 1..MAXPOL
!-----------------------------------------------------------------------------------------------
IF (NOPTST .EQ. 0 .AND. MAXPOL .GT. 0) THEN
FILE_POLY = TRIM (PATH_BIN) // 'POLYAD.BIN'
OPEN (UNIT = NSCR, FILE = FILE_POLY, DISPOSE = 'DELETE', &
& ACTION = 'READWRITE', FORM = 'UNFORMATTED')
IF (MAXQUA .EQ. 0) MAXQUA = MAXPOL
WRITE (NOUT,1410) MAXPOL,MAXQUA,(KPOLY(I),I=1,NQ)
WRITE (NOUT,1420)
NOPTST = 0
NSTATE = 0
MSIZE = 0 ! Maximum size of polyad blocks
DO 250 MODE = 0, 1
K = 0
DO 240 NPOLY = 1, MAXPOL
!---- Create states for given NPOLY and save to file NSCR
IOUT = 0 ! = MODE
CALL PSTATE (NQ,W,KPOLY,NPOLY,MAXQUA,MODE, & !ISYMSP,
& NSIZE,JQUANT,MAXTOT,LIMIT,IOUT)
IF (NSIZE .EQ. 0) CYCLE
IF (MODE .EQ. 0) THEN
MSIZE = MAX (MSIZE, NSIZE)
NSTATE = NSTATE + NSIZE
CYCLE
ENDIF
DO 230 I = 1, NSIZE
K = K + 1
DO Q = 1, NQ
IQUANT(Q,K) = JQUANT(Q,I)
END DO
!ISP = ISYMSTAT (NQ,IQUANT(1,K),ISYMSP)
!CALL VIBRENER (MDIM,NQ,IQUANT(1,K),W,XCON,EZERO,EHARM,EGROU,EVPT2)
CALL ENERHAOS (NQ,W,IQUANT(1,K),EZERO,EHARM)
WRITE (UNIT = QUANID, FMT = "('(',24(I2:','))") (IQUANT(Q,K),Q=1,NQ)
QUANID = TRIM(QUANID) // ')'
WRITE (NOUT,1430) K,NPOLY,I, & !LSYMSP(ISYMSP(ISP)),
& EHARM-EZERO,QUANID
230 CONTINUE
WRITE (NOUT,1440)
240 CONTINUE
IF (MODE .EQ. 0) THEN
ALLOCATE (IQUANT(NQ,0:NSTATE),JQUANT(NQ,0:MSIZE))
JQUANT = 0
IQUANT = 0 ! Ground state is initialized here
K = 0
I = 1
NPOLY = 0
!ISP = ISYMSTAT (NQ,IQUANT(1,0),ISYMSP)
!CALL VIBRENER (MDIM,NQ,IQUANT(1,K),W,XCON,EZERO,EHARM,EGROU,EVPT2)
CALL ENERHAOS (NQ,W,IQUANT(1,K),EZERO,EHARM)
WRITE (UNIT = QUANID, FMT = "('(',24(I2:','))") (IQUANT(Q,K),Q=1,NQ)
QUANID = TRIM(QUANID) // ')'
WRITE (NOUT,1430) K,NPOLY,I, & !LSYMSP(ISYMSP(ISP)),
& EHARM-EZERO,QUANID
WRITE (NOUT,1440)
ENDIF
IF (MODE .EQ. 1) DEALLOCATE (JQUANT)
250 CONTINUE
CLOSE (NSCR)
ENDIF
1410 FORMAT (/80('=')//'Generate states itemized by polyad number ', &
& '--'// &
& 'Maximum value of the polyad number =',I4/ &
& 'Maximum quanta of state excitation =',I4/ &
& 'Declared polyad coefficients: ',(15I4))
1420 FORMAT (/'Vibrational states itemized by polyad blocks and ', &
& 'symmetry species'/80('-')/ &
& ' ## Pol. Loc. Harm.Ener. ', &
& 'Vibrational State Quantum Numbers'/80('-'))
1430 FORMAT (I4,I5,I5,2X,F12.4,3X,A) !(I4,I5,I5,2X,A4,2F12.4,3X,A)
1440 FORMAT (80('-'))
!-----------------------------------------------------------------------------------------------
!=======================================================================
! Vibrational Configuration Interaction (VCI) Method
!=======================================================================
WRITE (NOUT,1500)
! ALLOCATE (ENERHARM(0:NSTATE),ENERVPT2(0:NSTATE), STAT = IERR) ! It will be used after VCI routine
! LOCERR = 6; IF (IERR .NE. 0) GO TO 990
! ALLOCATE (ENERHO(0:NSTATE),ENERCI(0:NSTATE),ENERPT(0:NSTATE), &
! & ENERPH(0:NSTATE), STAT = IERR)
! LOCERR = 6; IF (IERR .NE. 0) GO TO 990
ALLOCATE (NQUPOL(0:NSTATE), STAT = IERR)
LOCERR = 5; IF (IERR .NE. 0) GO TO 990
MAXPOL = 0
DO 300 L = 0, NSTATE
N = 0
DO Q = 1, NQ
N = N + KPOLY(Q) * IQUANT(Q,L)
END DO
NQUPOL(L) = N ! Polyad Quantum Number for each IQUANT state.
MAXPOL = MAX (MAXPOL, NQUPOL(L))
300 CONTINUE
IF ( NQEX + NQIN .NE. NQ) GOTO 950
WRITE (NOUT,1510) NQ,NQEX,NQIN,MPES,LPES(MPES),TOLCUT,NTHAM, & ! MODREP,
& NOPTST,NSTATE,NSTART,MAXPOL,MAXVRT,MULMOD,ENERMAX,MAX_PT, & ! KIND
& KIND(ONE) !,MANT
WRITE (* ,1510) NQ,NQEX,NQIN,MPES,LPES(MPES),TOLCUT,NTHAM, & ! MODREP,
& NOPTST,NSTATE,NSTART,MAXPOL,MAXVRT,MULMOD,ENERMAX,MAX_PT, & ! KIND
& KIND(ONE) !,MANT
IF (NSTART .GT. 0) THEN
WRITE (*,1520) NSTART,NSTATE
READ (*,"(I)") KEY
IF (KEY .EQ. 0) NSTART = 0
ENDIF
1500 FORMAT (/'Vibrational Rayleigh-Schrodinger Perturbation Theory:'/ &
& 'Calculate Ground State Level and Fundamental States.')
1510 FORMAT (/'RSPT parameters:'/ &
& 'The number of degrees of freedom = ',I8/ &
& 'The number of External Quant.Num.(v) = ',I8/ &
& 'The number of Internal Quant.Num.(l) = ',I8/ &
& 'The order of potential function = ',I8,2X,A16/ &
! & '(n)MR-PES (n-Mode Representation), n = ',I8/ &
& 'Tolerance for CVPT operator terms = ',ES8.1,' 1/Cm'/ &
& 'The number of operator terms in H = ',I8/ &
& 'Vibrational states set option, 1-4 = ',I8/ &
& 'The number of vibrational states = ',I8/ &
& 'Declared resumed calculation from = ',I8/ &
& 'Maximum polyad quantum number = ',I8/ &
& 'Max mode excitation of virtual states = ',I8/ &
& 'Number of modes for virtual states = ',I8/ &
& 'Energy limit for the virtual states = ',F10.1/ &
& 'Maximum Order of Perturbation Theory = ',I8/ &
& 'Numerical accuracy: KIND(REAL) = ',I8/ )
! & 'Numerical accuracy: Multi Precision = ',I8)
1520 FORMAT (/'Declaration: Use saved data about',I4,' states of',I4/ &
& 'Please acknowledge (1) or cancel (0) this option: '\)
IF (NQEX .NE. 0) DEALLOCATE(KPOLY)
!-----------------------------------------------------------------------------------------------
CALL LADDCOEF
KDIAG = 0 ! Eigenvalue Problems #
KSCAN = 0 ! Hamiltonian scans
KELEM = 0 ! Matrix element evals
IOPER = 0 ! Value[*,0] is not used
ALLOCATE (ISTATE(NQ))
!-----------------------------------------------------------------------------------------------
! Create virtual states: single, double, triple & quadruple
! (de-)excitations
!-----------------------------------------------------------------------------------------------
WRITE (NOUT,1600) MULMOD,MAXVRT,ENERMAX
IF (LOUT .GT. 0) WRITE (NOUT,1610)
DO 400 IPASS = 1, 2
IF (IPASS .EQ. 2) THEN
DO 310 V = 0, NSTATE
DO 310 Q = 1, NQ
IVIRTU(Q,V) = IQUANT(Q,V)
310 CONTINUE
ENDIF
NVIRT = NSTATE
!---- Generate virtual states: single (de-)excited configurations
LOCAL = 0
DO 340 NP = 1, NQ
DO K = 1, MAXVRT
!---- Create single (de-)excited configuration
ISTATE = 0
ISTATE(NP) = K
!---- Skip polyad states
DO V = 1, NSTATE
J = 0
DO I = 1, NQ
IF (ISTATE(I) .EQ. IQUANT(I,V)) J = J + 1
END DO
IF (J .EQ. NQ) EXIT
END DO
IF (J .EQ. NQ) CYCLE
CALL ENERHAOS (NQ,W,ISTATE,EZERO, ESTATE)
ESTATE = ESTATE - EZERO
IF (ESTATE .GT. ENERMAX) EXIT
NVIRT = NVIRT + 1
LOCAL = LOCAL + 1
IF (IPASS .EQ. 1) CYCLE
DO J = 1, NQ
IVIRTU(J,NVIRT) = ISTATE(J)
END DO
!---- Used for ID labeling only
NQUANT = IABS (K)
IF (LOUT .GT. 0 .AND. LOCAL .LE. LIMIT) THEN
WRITE (UNIT = QUANID, FMT = "('(',24(I2:','))") (ISTATE(J),J=1,NQ)
QUANID = TRIM(QUANID) // ')'
WRITE (NOUT,1620) NVIRT,SDTQ(1),NP,K,NQUANT,ESTATE,QUANID !
END IF
END DO
340 CONTINUE
NVIRT1 = LOCAL
1620 FORMAT (I4,3X,A6, 1X,I3,':',I2,22X,'=',I4,F13.2,4X,A80)
!---- Generate virtual states: double (de-)excited configurations
LOCAL = 0
DO 350 NP = 1, NQ - 1
DO 350 MQ = NP + 1, NQ
DO K = 1, MAXVRT
DO L = 1, MAXVRT
!---- Create double (de-)excited configuration
ISTATE = 0
ISTATE(NP) = K
ISTATE(MQ) = L
!---- Skip polyad states
DO V = 1, NSTATE
J = 0
DO I = 1, NQ
IF (ISTATE(I) .EQ. IQUANT(I,V)) J = J + 1
END DO
IF (J .EQ. NQ) EXIT
END DO
IF (J .EQ. NQ) CYCLE
CALL ENERHAOS (NQ,W,ISTATE,EZERO, ESTATE)
ESTATE = ESTATE - EZERO
IF (ESTATE .GT. ENERMAX) EXIT
NVIRT = NVIRT + 1
LOCAL = LOCAL + 1
IF (IPASS .EQ. 1) CYCLE
DO J = 1, NQ
IVIRTU(J,NVIRT) = ISTATE(J)
END DO
NQUANT = IABS (K) + IABS (L)
IF (LOUT .GT. 0 .AND. LOCAL .LE. LIMIT) THEN
WRITE (UNIT = QUANID, FMT = "('(',24(I2:','))") (ISTATE(J),J=1,NQ)
QUANID = TRIM(QUANID) // ')'
WRITE (NOUT,1630) NVIRT,SDTQ(2),NP,K,MQ,L,NQUANT,ESTATE,QUANID
END IF
END DO; END DO
350 CONTINUE
NVIRT2 = LOCAL
1630 FORMAT (I4,3X,A6,1X,I3,':',I2,',',I3,':',I2,15X,'=',I4,F13.2,4X,A80)
!---- Generate virtual states: triple (de-)excited configurations
NVIRT3 = 0
IF (MULMOD .LE. 2) GO TO 370
LOCAL = 0
DO 360 NP = 1, NQ - 2
DO 360 MQ = NP + 1, NQ - 1
DO 360 IR = MQ + 1, NQ
DO K = 1, MAXVRT
DO L = 1, MAXVRT
DO M = 1, MAXVRT
!---- Create triple (de-)excited configuration
ISTATE = 0
ISTATE(NP) = K
ISTATE(MQ) = L
ISTATE(IR) = M
!---- Skip polyad states
DO V = 1, NSTATE
J = 0
DO I = 1, NQ
IF (ISTATE(I) .EQ. IQUANT(I,V)) J = J + 1
END DO
IF (J .EQ. NQ) EXIT
END DO
IF (J .EQ. NQ) CYCLE
CALL ENERHAOS (NQ,W,ISTATE,EZERO, ESTATE)
ESTATE = ESTATE - EZERO
IF (ESTATE .GT. ENERMAX) EXIT
NVIRT = NVIRT + 1
LOCAL = LOCAL + 1
IF (IPASS .EQ. 1) CYCLE
DO J = 1, NQ
IVIRTU(J,NVIRT) = ISTATE(J)
END DO
NQUANT = IABS (K) + IABS (L) + IABS (M)
IF (LOUT .GT. 0 .AND. LOCAL .LE. LIMIT) THEN
WRITE (UNIT = QUANID, FMT = "('(',24(I2:','))") (ISTATE(J),J=1,NQ)
QUANID = TRIM(QUANID) // ')'
WRITE (NOUT,1640) NVIRT,SDTQ(3),NP,K,MQ,L,IR,M,NQUANT,ESTATE,QUANID
END IF
END DO; END DO; END DO
360 CONTINUE
NVIRT3 = LOCAL
1640 FORMAT (I4,3X,A6,1X,2(I3,':',I2,','),I3,':',I2,8X,'=',I4,F13.2,4X,A80)
!---- Generate virtual states: quadruple (de-)excited configurations
370 CONTINUE
NVIRT4 = 0
IF (MULMOD .LE. 3) GO TO 390
LOCAL = 0
DO 380 NP = 1, NQ - 3
DO 380 MQ = NP + 1, NQ - 2
DO 380 IR = MQ + 1, NQ - 1
DO 380 JS = IR + 1, NQ
DO K = 1, MAXVRT
DO L = 1, MAXVRT
DO M = 1, MAXVRT
DO N = 1, MAXVRT
!---- Create quadruple (de-)excited configuration
ISTATE = 0
ISTATE(NP) = K
ISTATE(MQ) = L
ISTATE(IR) = M
ISTATE(JS) = N
!---- Skip polyad states
DO V = 1, NSTATE
J = 0
DO I = 1, NQ
IF (ISTATE(I) .EQ. IQUANT(I,V)) J = J + 1
END DO
IF (J .EQ. NQ) EXIT
END DO
IF (J .EQ. NQ) CYCLE
CALL ENERHAOS (NQ,W,ISTATE,EZERO, ESTATE)
ESTATE = ESTATE - EZERO
IF (ESTATE .GT. ENERMAX) EXIT
NVIRT = NVIRT + 1
LOCAL = LOCAL + 1
IF (IPASS .EQ. 1) CYCLE
DO J = 1, NQ
IVIRTU(J,NVIRT) = ISTATE(J)
END DO
NQUANT = IABS (K) + IABS (L) + IABS (M) + IABS (N)
IF (LOUT .GT. 0 .AND. LOCAL .LE. LIMIT) THEN
WRITE (UNIT = QUANID, FMT = "('(',24(I2:','))") (ISTATE(J),J=1,NQ)
QUANID = TRIM(QUANID) // ')'
WRITE (NOUT,1650) NVIRT,SDTQ(3),NP,K,MQ,L,IR,M,JS,N,NQUANT,ESTATE,QUANID
END IF
END DO; END DO; END DO; END DO
380 CONTINUE
NVIRT4 = LOCAL
1650 FORMAT (I4,3X,A6,1X,3(I3,':',I2,','),I3,':',I2,' =',I4,F13.2,4X,A80)
390 IF (IPASS .EQ. 1) THEN
ALLOCATE (IVIRTU(NQ,0:NVIRT), STAT = IERR)
LOCERR = 7; IF (IERR .NE. 0) GO TO 990
! ALLOCATE (EVIRTU(0:NVIRT), STAT = IERR) ! It will be used in RSPT routine
! LOCERR = 8; IF (IERR .NE. 0) GO TO 990
ENDIF
400 CONTINUE ! IPASS = 1, 2
WRITE (NOUT,1660) NVIRT1,NVIRT2,NVIRT3,NVIRT4,NVIRT
1600 FORMAT (/'Generate virtual states for the Hamiltonian matrix:'/ & ! 2400
& 'Maximum number of different modes, MULMOD =',I8/ &
& 'Maximum mode quanta of excitation, MAXVRT =',I8/ &
& 'Energy limit for a virtual state, EnerMax =',F8.1/64('=')/ &
& 'Conditions: Energy of a virtual state <= EnerMax, and'/ &
& ' Max excitation on a single mode <= MAXVRT.')
1610 FORMAT (/'Virtual States Specifications:'/80('-')/ &
& ' No. Type Mode Excitations Total Energy', &
& ' Quantum Numbers'/80('-'))
1660 FORMAT (80('-')/ &
& 'Total number of 1-mode virtual states =',I8/ &
& 'Total number of 2-mode virtual states =',I8/ &
& 'Total number of 3-mode virtual states =',I8/ &
& 'Total number of 4-mode virtual states =',I8/48('=')/ &
& 'Total number of created virtual states =',I8)
!-----------------------------------------------------------------------------------------------
! Generate Lz Matrix. Note that Lz Matrix is Always Complex
! LZELE and ELEMAT is DOUBLE COMPLEX Scalar
! LZMAT is DOUBLE COMPLEX Matrix
!-----------------------------------------------------------------------------------------------
MVIRT = NVIRT + 1
ALLOCATE( IOPBRA(NQ), IOPKET(NQ) )
ALLOCATE( LZMAT (MVIRT,MVIRT), STAT = IERR)
LOCERR = 9; IF (IERR .NE. 0) GO TO 990
ICOUNT = 0
MCOUNT = MVIRT * (MVIRT + 1) / 2
CALL LWATCH (LSPLIT1)
ICENT0 = 0
CALL WATCH (ISPLIT)
WRITE (*,405) NVIRT,NTLZ
CALL COUNTER (0,0,1)
405 FORMAT (/'Evaluate L_z matrix:',I8,' states,',I8,' terms ...'\)
DO 420 K = 0, NVIRT
DO 420 L = 0, K
ICOUNT = ICOUNT + 1
ICENT1 = 100 * ICOUNT / MCOUNT
IF (ICENT1 .NE. ICENT0) THEN
CALL COUNTER (1,ICENT1,1)
ICENT0 = ICENT1
ENDIF
DO 410 Q = 1, NQ
IOPBRA(Q) = IVIRTU(Q,K) + 1
IOPKET(Q) = IVIRTU(Q,L) + 1
410 CONTINUE
ELEMAT = LZELE (NQ,NTLZ,LCOEF,LADLZ,IOPBRA,IOPKET)
LZMAT(K+1,L+1) = ELEMAT
LZMAT(L+1,K+1) = ELEMAT
420 CONTINUE
CALL COUNTER (1,-1,1)
CALL WATCH2 (ISPLIT)
DEALLOCATE (LCOEF, LADLZ)
IF (LOUT .EQ. 2) THEN
WRITE(NOUT,"(/' Lz Matrix has the following form'/)")
DO J = 1, MVIRT
WRITE (NOUT,"(200(1X,F10.5))") (IMAG(LZMAT(I,J)), I=1,MVIRT)
ENDDO
ENDIF
!--------------------------------------------
! LAPACK is Used in the VCI PROCEDURE
!
! ZHETRD - Reduces a Complex Hermitian Matrix to Tridiagonal Form.
! call Zhetrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
! ZUNGTR - Generate Unitary Matrix that Transfrom Hermitian Matrix to
! Tridiagonal Form
! call zungtr(uplo, n, a, lda, tau, work, lwork, info)
!
! ZSTEQR - Computes all eigenvalues and eigenvectors of a symmetric
! or Hermitian matrix reduced to tridiagonal form (QR algorithm).
! call zsteqr(compz, n, d, e, z, ldz, work, info)
! ZGEMM Computes a Matrix-Matrix product with Complex matrices.
!--------------------------------------------
ALLOCATE (TRIDIAG(MVIRT,MVIRT))
ALLOCATE (TAU(MVIRT-1) ,WORK(MVIRT) )
ALLOCATE (OFFDIAG(MVIRT) )
ALLOCATE (DIAG(MVIRT) )
ALLOCATE (EIGVECLZ(MVIRT,MVIRT), AUX(MVIRT,MVIRT))
CALL ZHETRD('U', MVIRT, LZMAT, MVIRT, DIAG, OFFDIAG, TAU, WORK, MVIRT, INFO)
AUX = LZMAT
CALL ZUNGTR('U', MVIRT, AUX, MVIRT, TAU, WORK, MVIRT, INFO)
CALL ZSTEQR('I', MVIRT, DIAG, OFFDIAG, LZMAT, MVIRT, WORK,INFO)
! EIGVECLZ = MATMUL(AUX, LZMAT)
CALL ZGEMM('N', 'N', MVIRT, MVIRT, MVIRT, CONE, AUX, MVIRT, LZMAT, MVIRT, CZERO, EIGVECLZ, MVIRT)
LZMAT = EIGVECLZ
DEALLOCATE (TRIDIAG, TAU, OFFDIAG, WORK, AUX, EIGVECLZ)
!-------------------------------------------------------------------------------------------
! MATRIX OUTPUT
!-------------------------------------------------------------------------------------------
IF ( LOUT .EQ. 2) THEN
WRITE(NOUT,"(/'Eigenvalues of Lz : '/)")
WRITE (NOUT, "(200(1X,F10.5))") (DIAG(I), I=1,MVIRT)
WRITE(NOUT,"(/'The Real Part of Lz Eigenvector : '/)")
DO J = 1, MVIRT
WRITE (NOUT,"(200(1X,F10.5))") (REAL(LZMAT(I,J)), I=1,MVIRT)
ENDDO
WRITE(NOUT,"(/'The Image Part of Lz Eigenvector : '/)")
DO J = 1, MVIRT
WRITE (NOUT,"(200(1X,F10.5))") (IMAG(LZMAT(I,J)), I=1,MVIRT)
ENDDO
ENDIF
!-------------------------------------------------------------------------------------------
! Select Eigenvectors Satisfying The Sayvetz Condition
! Use Eigenvectors to Form a Projector (PROJECT) and its Hermitian Conjugate (PROJTRAN),
! which will be used to transform Hamiltonian Matrix
!
! EIGVL Should be set to Chose Eigenvalue of Lz
!-------------------------------------------------------------------------------------------
WRITE(NOUT,*)
ALLOCATE(COPY1(MVIRT,1))
IF (NQEX .EQ. 0) THEN
DO MODE = 0, 1
IVIRTUL = 0
DO I = 1, MVIRT
IVIRTUL = IVIRTUL + 1 !
IF (MODE .EQ. 0 .AND. LOUT .EQ. 2) THEN !
WRITE (NOUT,"( F10.5)") DIAG(I) !
WRITE (NOUT,"(200(1X,F10.5))") (REAL(LZMAT(J,I)), J=1,MVIRT) !
WRITE (NOUT,"(200(1X,F10.5))") (IMAG(LZMAT(J,I)), J=1,MVIRT) !
ENDIF !
IF (MODE .EQ. 1) THEN !
DO J = 1, MVIRT !
COPY1(J,1) = LZMAT(J, I) !
PROJECT (IVIRTUL, J) = COPY1(J,1) !
PROJTRAN(J, IVIRTUL) = CONJG(COPY1(J,1)) !
ENDDO !
ENDIF ! !
ENDDO