-
Notifications
You must be signed in to change notification settings - Fork 0
/
ANHARM.FORT11
2572 lines (2572 loc) · 76.3 KB
/
ANHARM.FORT11
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
SUBROUTINE FENTRY(CC,IA,MAXCOR)
C***********************************************************************
C* MODIFICATION FOR FENTRY *
C* BY: YUKIO YAMAGUCHI *
C* DATE: FEBRUARY 16, 1989 *
C***********************************************************************
C* LAST UPDATED ON JANUARY 19, 1989 BY YUKIO YAMAGUCHI*
C***********************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
C2-16-89 DIMENSION CC(360000),IA(1)
DIMENSION CC(MAXCOR),IA(1)
DIMENSION I30(200),COORD(3,50)
DIMENSION I48(1024)
LOGICAL ISOMOL
CHARACTER*1 RAXIS
CHARACTER*25 RTYPE
DIMENSION ATM(36)
COMMON/VIB101/NATOM,N3N,NATRI,ILIN,NVIB
COMMON/VIB102/ITHREE,IFOUR,N3TOT,N4TOT
COMMON/VIB103/PARA,WAVE,CONST,CYCL,CONV
COMMON/VIB104/AIXX,AIYY,AIZZ,SUMW
COMMON/VIB105/FACT3,FACT4
COMMON/VIB106/ROTAA(3),ROTGC(3),ROTCM(3),ROTMH(3)
COMMON/VIB107/ITOP,IVSGL,IVDBL,IVTRP
COMMON/VIB108/IAXIS(3),NDEG(150),NDAB(150,5),IMAG(150)
COMMON/VIB109/CLIMIT,FLIM1,FLIM2
COMMON/VIB110/ISIGMA
COMMON/VIB201/CHG(50),X(50),Y(50),Z(50),W(50)
COMMON/VIB202/R(1275)
COMMON/VIB203/IOFF(150),IPRNT
COMMON/VIB204/SQM(150),ROOT(150),FREQ(150)
COMMON/VIB205/IFREQ,NFRQ(150)
COMMON/VIB206/RAXIS(3),RTYPE(6)
COMMON/VIB207/TABCD(3,3,3,3),PABC(3,3,3)
COMMON/VIB208/TAABB(3,3),TABAB(3,3)
COMMON/VIB209/SALP(3)
C2-16-89 EQUIVALENCE (CC,IA)
DATA ATM / 1.007825D+00 , 4.00260D+00 , 7.01600D+00 ,
* 9.01218D+00 , 11.00931D+00 , 12.00000D+00 ,
* 14.00307D+00 , 15.99491D+00 , 18.99840D+00 ,
* 19.99244D+00 , 22.98977D+00 , 23.98504D+00 ,
* 26.98154D+00 , 27.97693D+00 , 30.97376D+00 ,
* 31.97207D+00 , 34.96885D+00 , 39.96238D+00 ,
* 38.96371D+00 , 39.96259D+00 , 44.95591D+00 ,
* 47.94795D+00 , 50.94396D+00 , 51.94051D+00 ,
* 54.93805D+00 , 55.93494D+00 , 58.93320D+00 ,
* 57.93535D+00 , 62.92960D+00 , 63.92915D+00 ,
* 68.92558D+00 , 73.92118D+00 , 74.92160D+00 ,
* 79.91652D+00 , 78.91834D+00 , 83.91151D+00 /
DATA PH,BK,AVN / 6.626176D+00 , 1.380662D+00 , 6.022045D+00 /
DATA A0,HE,EV / 0.52917706D+00 , 4.359813653D+00 , 1.6021892D+00 /
DATA CL,TEMP / 2.99792458D+00 , 300.0D+00 /
DATA PI,ALN / 3.1415926536D+00 , 2.3025850930D+00 /
1 FORMAT(//,2X,' THE ANHARMONIC CONSTANTS PROGRAM'/)
2 FORMAT(10I5)
3 FORMAT(F20.10)
4 FORMAT(//,2X,' PARAMETERS FOR THE CALCULATION'//
* 2X,' ISOTOP = ',I5/
* 2X,' IGEOMT = ',I5/
* 2X,' ITHREE = ',I5/
* 2X,' IFOUR = ',I5/
* 2X,' IFREQ = ',I5/
* 2X,' ICORIO = ',I5/
* 2X,' IFERM1 = ',I5/
* 2X,' IFERM2 = ',I5/
* 2X,' ISIGMA = ',I5/
* 2X,' IPRNT = ',I5/
* 2X,' NATOM = ',I5/
* 2X,' N3N = ',I5/
* 2X,' NATRI = ',I5/
* 2X,' N3TOT = ',I5/
* 2X,' N4TOT = ',I5/)
5 FORMAT(//,2X,' SCF ENERGY = ',F20.10/)
6 FORMAT(//,2X,' MOLECULAR GEOMETRY (IN BOHR)'/
1 5X,4H NO.,11X,2H X,18X,2H Y,18X,2H Z,18X,2H W/)
7 FORMAT(4F20.10)
8 FORMAT(2X,I7,5X,4F20.10)
9 FORMAT(//,2X,' MOLECULAR GEOMETRY (IN A)'/
1 5X,4H NO.,11X,2H X,18X,2H Y,18X,2H Z,18X,2H W/)
10 FORMAT(F20.10)
11 FORMAT(//,2X,' ICMAX = ',I10,5X,' MAXCOR = ',I10/)
12 FORMAT('1',//,2X,' THE ANHARMONIC CONSTANTS PROGRAM'/)
C
CALL TSTART(6)
CALL NOUNFL
C
C PHYSICAL CONSTANTS
PARA=1.0D+00/AVN
WAVE=1.0D+04*DSQRT(AVN)/(2.0D+00*PI*CL)
CONST=1.0D+02*(PH*AVN)/(8.0D+00*PI*PI*CL)
CYCL=1.0D+06*(PH*AVN)/(8.0D+00*PI*PI)
CONV=1.0D-01*(PH*CL)/(2.0D+00*BK*TEMP)
C ALAM IS IN AMU-1.A-2
ALAM=4.0D-02*(PI*PI*CL)/(PH*AVN)
C SQLAM IS IN AMU-1/2.A-1
SQLAM=DSQRT(ALAM)
C FACT3 IS IN CM-1
FACT3=1.0D+06/(SQLAM*SQLAM*SQLAM*PH*CL)
C FACT4 IS IN CM-1
FACT4=1.0D+06/(ALAM*ALAM*PH*CL)
C
INPUT=5
ITAPE6=6
ITAP15=15
ITAP30=30
ITAP48=48
NSECT=1024
C2-16-89 MAXCOR=360000
C
CALL LOCATE(INPUT,'# ANHARM #',IERR)
C
WRITE(6,1)
C
C***INPUT PARAMETERS FOR THE CALCULATION***
C ISOTOP IS A PARAMETER FOR ISOTOPIC MOLECULES
C ISOTOP =0 AND 1 :USE REGULAR ATOMIC MASS
C ISOTOP =2 OR MORE :NUMBER OF ISOTOPIC MOLECULES
C IGEOMT IS A PARAMETER FOR INPUT DATA
C IGEOMT =0 :READ IN GEOMETRY FROM TAPE 30
C IGEOMT =5 :READ IN GEOMETRY FROM TAPE 5
C IGEOMT =48:READ IN GEOMETRY FROM FILE48
C ITHREE IS A PARAMETER FOR THIRD DERIVATIVES
C ITHREE =0 :DO NOT READ IN THIRD DERIVATIVES
C ITHREE =1 OR 20:READ IN THIRD DERIVATIVES FROM FILE20
C ITHREE =48:READ IN THIRD DERIVATIVES FROM FILE48
C IFOUR IS A PARAMETER FOR FOURTH DERIVATIVES
C IFOUR =0 :DO NOT READ IN FOURTH DERIVATIVES
C IFOUR =1 OR 24:READ IN FOURTH DERIVATIVES FROM FILE24
C IFOUR =48:READ IN FOURTH DERIVATIVES FROM FILE48
C IFREQ =0 :CALCULATE 3N-6(OR 3N-5) SETS OF ANHARMONIC CONSTANTS
C IFREQ =N :CALCULATE IFREQ SETS OF ANHARMONIC CONSTANTS
C ICORIO =0 :USE A DEFAULT VALU FOR CORIOLIS RESONANCE
C ICORIO =1 :READ IN A THRESHOLD VALU FOR CORIOLIS RESONANCE
C IFERM1 =0 :USE A DEFAULT VALU FOR TYPE 1 FERMI RESONANCE
C TYPE1----W(I)=W(J)
C IFERM1 =1 :READ IN A THRESHOLD VALU FOR TYPE 1 FERMI RESONANCE
C IFERM2 =0 :USE A DEFAULT VALU FOR TYPE 2 FERMI RESONANCE
C TYPE2---W(I)+W(J)=W(K)
C IFERM2 =1 :READ IN A THRESHOLD VALU FOR TYPE 2 FERMI RESONANCE
C ISIGMA =0 :USE A(0') VALUE TO EVALUATE SIGMA VALUE THAT IS NEEDED
C FOR CENTRIFUGAL DISTORTION CONSTANTS
C ISIGMA =1 :USE A(E) VALUE TO EVALUATE SIGMA VALUE THAT IS NEEDED
C FOR CENTRIFUGAL DISTORTION CONSTANTS
READ(5,2) ISOTOP,IGEOMT,ITHREE,IFOUR,IFREQ,ICORIO,IFERM1,IFERM2,
1 ISIGMA,IPRNT
IF(IGEOMT.EQ.0) IGEOMT=30
IF(ITHREE.EQ.1) ITHREE=20
IF(IFOUR.EQ.1) IFOUR=24
C
REWIND ITAP15
READ(ITAP15,2) NATOM,NSTORE
C
IOFF(1)=0
DO 101 I=1,149
IOFF(I+1)=IOFF(I)+I
101 CONTINUE
C
IF(IGEOMT.EQ.30) THEN
CALL RFILE(ITAP30)
CALL WREADW(ITAP30,I30,200,101,JUNK)
IEND=I30(1)
MPOINT=I30(2)
MCONST=I30(3)
MCALCS=I30(4)
NCALCS=I30(5)
NATOM=I30(19)
N3N=NATOM*3
JUNK=101+MCONST
CALL WREADW(ITAP30,IPT,1,JUNK,JUNK)
CALL WREADW(ITAP30,CHG,NATOM*2,IPT,JUNK)
IGEOP=100+MCONST+MPOINT+NCALCS
CALL WREADW(ITAP30,LOCCAL,1,IGEOP,JUNK)
IGEOP=LOCCAL+60
CALL WREADW(ITAP30,I30,20,IGEOP,IGEOP)
CALL WREADW(ITAP30,COORD,N3N*2,IGEOP,JUNK)
JUNK=JUNK+INTOWP(1)
CALL WREADW(ITAP30,ESCF,2,JUNK,JUNK)
CALL RCLOSE(ITAP30,3)
END IF
C
IF(IGEOMT.EQ.48) THEN
CALL RFILE(ITAP48)
CALL RREAD(ITAP48,I48,NSECT,2)
NATOM=I48(4)
N3N=I48(5)
CALL RREAD(ITAP48,COORD,N3N*2,5)
CALL RREAD(ITAP48,W,NATOM*2,6)
CALL RCLOSE(ITAP48,3)
END IF
C
N3N=NATOM*3
NATRI=IOFF(N3N+1)
N3TOT=N3N*(N3N+1)*(N3N+2)/6
N4TOT=N3N*(N3N+1)*(N3N+2)*(N3N+3)/24
C
C SET AN ARRAY FOR ANHARMONICITY CALCULATION
DO 102 I=1,N3N
102 NFRQ(I)=I
CALL LOCATE(INPUT,'# ANHARM #',IERR)
READ(5,2) ISOTOP,IGEOMT,ITHXXX,IFXXX,IFREQ,ICORIO,IFERM1,IFERM2,
1 ISIGMA,IPRNT
IF(ISOTOP.EQ.0) ISOTOP=1
IF(IFREQ.NE.0) THEN
NFREQ=IABS(IFREQ)
READ(5,2) (NFRQ(I),I=1,NFREQ)
END IF
IF(ICORIO.NE.0) THEN
READ(5,3) CLIMIT
ELSE
CLIMIT=100.0D+00
END IF
IF(IFERM1.NE.0) THEN
READ(5,3) FLIM1
ELSE
FLIM1=100.0D+00
END IF
IF(IFERM2.NE.0) THEN
READ(5,3) FLIM2
ELSE
FLIM2=100.0D+00
END IF
C
WRITE(6,4) ISOTOP,IGEOMT,ITHREE,IFOUR,IFREQ,ICORIO,IFERM1,IFERM2,
1 ISIGMA,IPRNT,NATOM,N3N,NATRI,N3TOT,N4TOT
C
WRITE(6,5) ESCF
WRITE(6,6)
DO 103 I=1,NATOM
IF(IGEOMT.EQ.5) GO TO 201
X(I)=COORD(1,I)
Y(I)=COORD(2,I)
Z(I)=COORD(3,I)
II=CHG(I)
IF(IGEOMT.EQ.48) GO TO 202
W(I)=ATM(II)
GO TO 202
201 CONTINUE
READ(5,7) COORD(1,I),COORD(2,I),COORD(3,I),W(I)
X(I)=COORD(1,I)
Y(I)=COORD(2,I)
Z(I)=COORD(3,I)
202 CONTINUE
WRITE(6,8) I,X(I),Y(I),Z(I),W(I)
103 CONTINUE
CALL DIST(1)
WRITE(6,9)
DO 104 I=1,NATOM
X(I)=X(I)*A0
Y(I)=Y(I)*A0
Z(I)=Z(I)*A0
WRITE(6,8) I,X(I),Y(I),Z(I),W(I)
104 CONTINUE
CALL DIST(1)
C
C DYNAMIC ALLOCATION
C IC1 : FX (N3N*N3N)
C IC2 : ELX (N3N*N3N)
C IC3 : ELXM (N3N*N3N)
C IC4 : UX (N3N*N3N)
C IC5 : EPX (N3N*N3N)
C IC6 : EPXM (N3N*N3N)
C IC7 : ZETA (N3N*N3N*3)
C IC8 : AXY (3*3*N3N)
C IC9 : BXY (3*3*N3N)
C IC10 : CXY (3*3*N3N)
C IC11 : BXYZ (3*3*3*N3N)
C IC12 : CXYZ (3*3*3*N3N)
C IC13 : F33 (N3TOT)
C IC14 : F3X (N3N*N3N*N3N)
C IC15 : F3Q (N3N*N3N*N3N)
C IC16 : F44 (N4TOT)
C IC17 : F4X (N3N*N3N*N3N*N3N)
C IC18 : F4Q (N3N*N3N*N3N*N3N)
C
C DYNAMIC ALLOCATION OF CORE MEMORY AND CLEARANCE OF IT
IC1=1
IC2=IC1+N3N*N3N
IC3=IC2+N3N*N3N
IC4=IC3+N3N*N3N
IC5=IC4+N3N*N3N
IC6=IC5+N3N*N3N
IC7=IC6+N3N*N3N
IC8=IC7+N3N*N3N*3
IC9=IC8+3*3*N3N
IC10=IC9+3*3*N3N
IC11=IC10+3*3*N3N
IC12=IC11+3*3*3*N3N
IC13=IC12+3*3*3*N3N
IC14=IC13+N3TOT
IC15=IC14+N3N*N3N*N3N
IC16=IC15+N3N*N3N*N3N
IC17=IC16+N4TOT
IC18=IC17+N3N*N3N*N3N*N3N
ICKEEP=IC18+N3N*N3N*N3N*N3N
CALL ZERO(CC(1),ICKEEP)
C
C FORM THE FORCE CONSTANT MATRIX
WRITE(3,21)
21 FORMAT(//,2X,' NOW YOU ARE IN FCONST'/)
ICMAX=ICKEEP
WRITE(3,11) ICMAX,MAXCOR
C.................FX......
CALL FCONST(CC(IC1))
C
C CALCULATE THE CENTER OF MASS AND MOMENTS OF INERTIA
C FOR THE PARENT MOLECULE
WRITE(3,22)
22 FORMAT(//,2X,' NOW YOU ARE IN MOMENT'/)
IC19=ICKEEP
ICMAX=IC19+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C.................UX EE.......
CALL MOMENT(CC(IC4),CC(IC19))
C
ISOMOL=.FALSE.
NISO=0
400 NISO=NISO+1
IF(NISO.EQ.1) GO TO 401
IF(NISO.GT.ISOTOP) GO TO 220
WRITE(6,12)
ISOMOL=.TRUE.
C
C READ IN ISOTOPIC ATOMIC MASS
WRITE(6,6)
DO 105 I=1,NATOM
READ(5,10) W(I)
X(I)=COORD(1,I)
Y(I)=COORD(2,I)
Z(I)=COORD(3,I)
WRITE(6,8) I,X(I),Y(I),Z(I),W(I)
105 CONTINUE
CALL DIST(1)
WRITE(6,9)
DO 106 I=1,NATOM
X(I)=X(I)*A0
Y(I)=Y(I)*A0
Z(I)=Z(I)*A0
WRITE(6,8) I,X(I),Y(I),Z(I),W(I)
106 CONTINUE
CALL DIST(1)
C
C CALCULATE THE CENTER OF MASS AND MOMENTS OF INERTIA
C FOR AN ISOTOPIC MOLECULE
CALL ZERO(CC(IC2),ICKEEP)
WRITE(3,22)
IC19=ICKEEP
ICMAX=IC19+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C.................UX EE.......
CALL MOMENT(CC(IC4),CC(IC19))
C
C THE NORMAL COORDINATE ANALYSIS
C
C THE FX MATRIX METHOD
401 CONTINUE
WRITE(3,23)
23 FORMAT(//,2X,' NOW YOU ARE IN VIBFX'/)
IC19=ICKEEP
IC20=IC19+N3N*N3N
IC21=IC20+N3N
IC22=IC21+N3N
IC23=IC22+N3N
IC24=IC23+NATRI
ICMAX=IC24+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C................FX ELX ELXM UX EPX EPXM....
CALL VIBFX(CC(IC1),CC(IC2),CC(IC3),CC(IC4),CC(IC5),CC(IC6),
C................EE VALU FV1 FV2 FXM DD.......
1 CC(IC19),CC(IC20),CC(IC21),CC(IC22),CC(IC23),CC(IC24),
C................CC.......
2 CC(IC24))
C
C CORIOLIS COUPLING CONSTANTS
WRITE(3,24)
24 FORMAT(//,2X,' NOW YOU ARE IN CORIOL'/)
IC19=ICKEEP
IC20=IC19+N3N*N3N
IC21=IC20+N3N*N3N
ICMAX=IC21+N3N*N3N*3
WRITE(3,11) ICMAX,MAXCOR
C.................EPXM ZETA EE FF AM.......
CALL CORIOL(CC(IC6),CC(IC7),CC(IC19),CC(IC20),CC(IC21))
C
C CALCULATE AXY MATRICES
WRITE(3,25)
25 FORMAT(//,2X,' NOW YOU ARE IN ABCMAT'/)
IC19=ICKEEP
ICMAX=IC19+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C.................EPXM AXY BXY CXY EE.......
CALL ABCMAT(CC(IC6),CC(IC8),CC(IC9),CC(IC10),CC(IC19))
C
C CALCULATE TAU MATRICES
WRITE(3,27)
27 FORMAT(//,2X,' NOW YOU ARE IN TAUMAT'/)
IC19=ICKEEP
ICMAX=IC19+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C.................AXY BXY.....
CALL TAUMAT(CC(IC8),CC(IC9))
C
C CALCULATE BXYZ MATRICES
WRITE(3,28)
28 FORMAT(//,2X,' NOW YOU ARE IN BCXYZ'/)
ICMAX=ICKEEP
WRITE(3,11) ICMAX,MAXCOR
C................ZETA BXY BXYZ CXYZ.....
CALL BCXYZ(CC(IC7),CC(IC9),CC(IC11),CC(IC12))
C
C READ CUBIC AND QUARTIC FORCE CONSTANTS FROM THE TAPE48
IF(ITHREE.NE.48.AND.IFOUR.NE.48) GO TO 402
WRITE(3,29)
29 FORMAT(//,2X,' NOW YOU ARE IN READ48'/)
IC19=ICKEEP
ICMAX=IC19+NVIB*NVIB*NVIB*NVIB
WRITE(3,11) ICMAX,MAXCOR
C.................F3Q F4Q F3X F4X DD CC.......
CALL READ48(CC(IC15),CC(IC18),CC(IC19),CC(IC19),CC(IC19),CC(IC19))
C
C CALCULATE CUBIC FORCE CONSTANTS
402 CONTINUE
IF(ITHREE.EQ.0) GO TO 303
IF(ITHREE.EQ.48) GO TO 303
WRITE(3,30)
30 FORMAT(//,2X,' NOW YOU ARE IN FORCE3'/)
IC19=ICKEEP
IC20=IC19+N3N*N3N
IC21=IC20+N3N*N3N
ICMAX=IC21+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C.................ELX F33 F3X F3Q EE DD.......
CALL FORCE3(CC(IC2),CC(IC13),CC(IC14),CC(IC15),CC(IC19),CC(IC20),
C.................GG.......
1 CC(IC21))
C
C CALCULATE PHI MATRICES
303 CONTINUE
WRITE(3,31)
31 FORMAT(//,2X,' NOW YOU ARE IN PHIMAT'/)
IC19=ICKEEP
ICMAX=IC19+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C.................BXY CXY BXYZ CXYZ F3Q EE.......
CALL PHIMAT(CC(IC9),CC(IC10),CC(IC11),CC(IC12),CC(IC15),CC(IC19))
C
C CALCULATE QUARTIC FORCE CONSTANTS
IF(IFOUR.EQ.0) GO TO 304
IF(IFOUR.EQ.48) GO TO 304
WRITE(3,32)
32 FORMAT(//,2X,' NOW YOU ARE IN FORCE4'/)
IC19=ICKEEP
IC20=IC19+N3N*N3N
IC21=IC20+N3N*N3N
ICMAX=IC21+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C.................ELX F44 F4X F4Q EE DD.......
CALL FORCE4(CC(IC2),CC(IC16),CC(IC17),CC(IC18),CC(IC19),CC(IC20),
C.................GG.......
1 CC(IC21))
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%ANHARMONICITY ANALYSIS %%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C SELECT THE ROTATIONAL TYPE OF THE MOLECULE
304 CONTINUE
GO TO (501,502,503,504,504,505),ITOP
C
C:::::::::::::::::::::::
C:::DIATOMIC MOLECULE:::
C:::::::::::::::::::::::
501 CONTINUE
WRITE(3,33)
33 FORMAT(//,2X,' NOW YOU ARE IN DIATOM'/)
C% IC19=ICKEEP
C% ICMAX=IC19+N3N*N3N
C% WRITE(3,11) ICMAX,MAXCOR
C%................F3Q F4Q EE.......
C% CALL DIATOM(CC(IC15),CC(IC18),CC(IC19))
WRITE(6,*) ' '
WRITE(6,*) ' ***************************************************'
WRITE(6,*) ' SUBROUTINE FOR A DIATOMIC MOLECULE IS NOT AVAILABLE'
WRITE(6,*) ' ***************************************************'
WRITE(6,*) ' '
GO TO 215
C
C::::::::::::::::::::::::::::::::
C:::LINEAR POLYATOMIC MOLECULE:::
C::::::::::::::::::::::::::::::::
502 CONTINUE
WRITE(3,34)
34 FORMAT(//,2X,' NOW YOU ARE IN LINEAR'/)
C% IC19=ICKEEP
C% IC20=IC19+N3N*N3N
C% ICMAX=IC20+N3N*N3N
C% WRITE(3,11) ICMAX,MAXCOR
C%................ZETA CXY F3Q CLTYP....
C% CALL LTYPLI(CC(IC7),CC(IC10),CC(IC15),CC(IC19))
C%................ZETA AXY F3Q EE RVCNST...
C% CALL ROVLIR(CC(IC7),CC(IC8),CC(IC15),CC(IC19),CC(IC20))
C%................BXY CXY F3Q EE.......
C% CALL CENTLI(CC(IC9),CC(IC10),CC(IC15),CC(IC19))
C% IF(IVDBL.LE.1) GO TO 205
C%................F3Q F4Q CLTYP....
C% CALL VLTYLI(CC(IC15),CC(IC18),CC(IC19))
C%205 CONTINUE
C%................ZETA F3Q F4Q XIJ GIJ......
C% CALL LINEAR(CC(IC7),CC(IC15),CC(IC18),CC(IC19),CC(IC20))
WRITE(6,*) ' '
WRITE(6,*) ' *************************************************'
WRITE(6,*) ' SUBROUTINE FOR A LINEAR MOLECULE IS NOT AVAILABLE'
WRITE(6,*) ' *************************************************'
WRITE(6,*) ' '
GO TO 215
C
C:::::::::::::::::::::::::::::
C:::ASYMMETRIC TOP MOLECULE:::
C:::::::::::::::::::::::::::::
503 CONTINUE
WRITE(3,35)
35 FORMAT(//,2X,' NOW YOU ARE IN ASYTOP'/)
IC19=ICKEEP
IC20=IC19+N3N*N3N
ICMAX=IC20+N3N*N3N
WRITE(3,11) ICMAX,MAXCOR
C.................ZETA AXY F3Q EE RVCNST...
CALL ROVASY(CC(IC7),CC(IC8),CC(IC15),CC(IC19),CC(IC20))
C.................F3Q EE.......
CALL CENTAS(CC(IC15),CC(IC19))
C.................ZETA F3Q F4Q XIJ......
CALL ASYTOP(CC(IC7),CC(IC15),CC(IC18),CC(IC19))
GO TO 215
C
C::::::::::::::::::::::::::::
C:::SYMMETRIC TOP MOLECULE:::
C::::::::::::::::::::::::::::
504 CONTINUE
WRITE(3,36)
36 FORMAT(//,2X,' NOW YOU ARE IN SYMTOP'/)
C% IC19=ICKEEP
C% IC20=IC19+N3N*N3N
C% ICMAX=IC20+N3N*N3N
C% WRITE(3,11) ICMAX,MAXCOR
C%................ZETA AXY F3Q EE CLTYP....
C% CALL LTYPSY(CC(IC7),CC(IC8),CC(IC15),CC(IC19),CC(IC20))
C%................ZETA AXY F3Q EE RVCNST...
C% CALL ROVSYM(CC(IC7),CC(IC8),CC(IC15),CC(IC19),CC(IC20))
C%................F3Q EE.......
C% CALL CENTSY(CC(IC15),CC(IC19))
C% IF(IVDBL.LE.1) GO TO 206
C%................ZETA AXY F3Q F4Q CLTYP....
C% CALL VLTYSY(CC(IC7),CC(IC8),CC(IC15),CC(IC18),CC(IC19))
C%206 CONTINUE
C%................ZETA F3Q F4Q XIJ GIJ......
C% CALL SYMTOP(CC(IC7),CC(IC15),CC(IC18),CC(IC19),CC(IC20))
WRITE(6,*) ' '
WRITE(6,*) ' ***********************************************'
WRITE(6,*) ' SUBROUTINE FOR A SYMMETRIC TOP IS NOT AVAILABLE'
WRITE(6,*) ' ***********************************************'
WRITE(6,*) ' '
GO TO 215
C
C::::::::::::::::::::::::::::
C:::SPHERICAL TOP MOLECULE:::
C::::::::::::::::::::::::::::
505 CONTINUE
WRITE(3,37)
37 FORMAT(//,2X,' NOW YOU ARE IN SPHTOP'/)
C% IC19=ICKEEP
C% IC20=IC19+N3N*N3N
C% ICMAX=IC20+N3N*N3N
C% WRITE(3,11) ICMAX,MAXCOR
C%................ZETA AXY F3Q EE RVCNST...
C% CALL ROVSPH(CC(IC7),CC(IC8),CC(IC15),CC(IC19),CC(IC20))
C%................F3Q EE.......
C% CALL CENTSP(CC(IC15),CC(IC19))
C%................ZETA F3Q F4Q XIJ GIJ......
C% CALL SPHTOP(CC(IC7),CC(IC15),CC(IC18),CC(IC19),CC(IC20))
WRITE(6,*) ' '
WRITE(6,*) ' ***********************************************'
WRITE(6,*) ' SUBROUTINE FOR A SPHERICAL TOP IS NOT AVAILABLE'
WRITE(6,*) ' ***********************************************'
WRITE(6,*) ' '
C
C##########################################
C RETURN FOR ANOTHER ISOTOPIC MOLECULE
C##########################################
215 CONTINUE
GO TO 400
C
220 CONTINUE
CALL TSTOP(6)
C
RETURN
C2-16-89 STOP
END
SUBROUTINE DIST(IPRNT)
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/VIB101/NATOM,N3N,NATRI,ILIN,NVIB
COMMON/VIB201/CHG(50),X(50),Y(50),Z(50),W(50)
COMMON/VIB202/R(1275)
DATA A00 / 0.0D+00 /
1 FORMAT(//,2X,' INTERATOMIC DISTANCE MATRIX'/)
C
IJ=0
DO 101 I=1,NATOM
DO 101 J=1,I
IJ=IJ+1
R(IJ)=A00
IF(I.EQ.J) GO TO 101
XD=X(I)-X(J)
YD=Y(I)-Y(J)
ZD=Z(I)-Z(J)
R(IJ)=DSQRT(XD*XD+YD*YD+ZD*ZD)
101 CONTINUE
IF(IPRNT.EQ.0) GO TO 201
WRITE(6,1)
CALL PRINT(R,1275,NATOM,6)
201 RETURN
END
SUBROUTINE FCONST(FX)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION FX(N3N,N3N)
COMMON/VIB101/NATOM,N3N,NATRI,ILIN,NVIB
DATA HE,A0 / 4.359813653D+00 , 0.52917706D+00 /
DATA ITAP15 / 15 /
1 FORMAT(2I5)
2 FORMAT(3F20.10)
3 FORMAT(//,2X,' FX MATRIX (IN HARTREE/BOHR+2)'/)
4 FORMAT(//,2X,' FX MATRIX (IN MDYN/A)'/)
C
C FUNIT IS IN MDYN/A
FUNIT=HE/(A0*A0)
C
C READ IN ANALYTICAL SECOND DERIVATIVES
C WITH RESPECT TO NUCLEAR COORDINATE
REWIND ITAP15
READ(ITAP15,1) NNATOM,NNSTOR
READ(ITAP15,2) ((FX(I,J),J=1,N3N),I=1,N3N)
REWIND ITAP15
WRITE(6,3)
CALL MATOUT(FX,N3N,N3N,N3N,N3N,6)
DO 101 I=1,N3N
DO 101 J=1,N3N
FX(I,J)=FX(I,J)*FUNIT
101 CONTINUE
WRITE(6,4)
CALL MATOUT(FX,N3N,N3N,N3N,N3N,6)
C
RETURN
END
SUBROUTINE MOMENT(UX,EE)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION UX(N3N,N3N),EE(N3N,N3N)
DIMENSION T(6),PT(3),EG(3,3),FV1(3),FV2(3)
CHARACTER*1 RAXIS
CHARACTER*25 RTYPE
COMMON/VIB101/NATOM,N3N,NATRI,ILIN,NVIB
COMMON/VIB103/PARA,WAVE,CONST,CYCL,CONV
COMMON/VIB104/AIXX,AIYY,AIZZ,SUMW
COMMON/VIB106/ROTAA(3),ROTGC(3),ROTCM(3),ROTMH(3)
COMMON/VIB107/ITOP,IVSGL,IVDBL,IVTRP
COMMON/VIB108/IAXIS(3),NDEG(150),NDAB(150,5),IMAG(150)
COMMON/VIB201/CHG(50),X(50),Y(50),Z(50),W(50)
COMMON/VIB206/RAXIS(3),RTYPE(6)
DATA A00,ONE,TWO / 0.0D+00 , 1.0D+00 , 2.0D+00 /
DATA A0 / 0.52917706D+00 /
DATA DLIMIT / 1.0D-02 /
DATA RLIMIT / 1.0D-05 /
1 FORMAT(//,2X,' CENTER OF MASS'/2X,' CMX',11X,' CMY',11X,' CMZ'//
1 2X,3(5X,F10.7)/)
2 FORMAT(//,2X,' CARTESIAN COORDINATES W.R.T. CENTER OF MASS (IN A)'
1/5X,4H NO.,11X,2H X,18X,2H Y,18X,2H Z,18X,2H W/)
3 FORMAT(2X,I7,5X,4F20.10)
4 FORMAT(//,2X,' INERTIA TENSOR'/)
5 FORMAT(//,2X,' PRINCIPAL MOMENTS OF INERTIA AT EQUILIBRIUM GEOMETR
1Y'/
2 4X,5H AXIS,9X,11H IN AMU.A+2,9X,10H IN G.CM+2,10X,
3 8H IN CM-1,15X,7H IN MHZ/38X,8H (*D+39)/)
6 FORMAT(8X,A1,3X,3F20.10,F20.5)
7 FORMAT(//,2X,' EIGENVECTOR OF INERTIA TENSOR'/)
8 FORMAT(//,2X,' THE ROTATIONAL TYPE OF THIS ISOTOPIC MOLECULE = ',
1 2X,A25/)
9 FORMAT(/,2X,' ASYMMETRY PARAMETERS'/
1 2X,' B(P) = ',F10.5/
2 2X,' B(O) = ',F10.5/
3 2X,' KAPPA = ',F10.5/)
10 FORMAT(//,2X,' UX MATRIX'/)
11 FORMAT(//,2X,' UX * UX(T) = E '/)
12 FORMAT(//,2X,' CARTESIAN COORDINATES W.R.T. PRINCIPAL AXIS (IN BOH
1R)'/5X,4H NO.,11X,2H X,18X,2H Y,18X,2H Z,18X,2H W/)
13 FORMAT(//,2X,' CARTESIAN COORDINATES W.R.T. PRINCIPAL AXIS (IN A)'
1/5X,4H NO.,11X,2H X,18X,2H Y,18X,2H Z,18X,2H W/)
C
C THE CALCULATION OF CENTER OF MASS
SUMW=A00
SUMWX=A00
SUMWY=A00
SUMWZ=A00
DO 101 I=1,NATOM
SUMW=SUMW+W(I)
SUMWX=SUMWX+W(I)*X(I)
SUMWY=SUMWY+W(I)*Y(I)
SUMWZ=SUMWZ+W(I)*Z(I)
101 CONTINUE
CMX=SUMWX/SUMW
CMY=SUMWY/SUMW
CMZ=SUMWZ/SUMW
WRITE(6,1) CMX,CMY,CMZ
C THE CARTESIAN COORDINATES W.R.T. CENTER OF MASS
WRITE(6,2)
DO 102 I=1,NATOM
X(I)=X(I)-CMX
Y(I)=Y(I)-CMY
Z(I)=Z(I)-CMZ
WRITE(6,3) I,X(I),Y(I),Z(I),W(I)
102 CONTINUE
C
C THE CALCULATION OF INERTIA TENSOR
CALL ZERO(T,6)
DO 104 I=1,NATOM
T(1)=T(1)+W(I)*(Y(I)**2+Z(I)**2)
T(3)=T(3)+W(I)*(Z(I)**2+X(I)**2)
T(6)=T(6)+W(I)*(X(I)**2+Y(I)**2)
T(2)=T(2)-W(I)*X(I)*Y(I)
T(4)=T(4)-W(I)*Z(I)*X(I)
T(5)=T(5)-W(I)*Y(I)*Z(I)
104 CONTINUE
AIXX=T(1)
AIYY=T(3)
AIZZ=T(6)
ILIN=0
IF(AIXX.LE.DLIMIT) ILIN=1
IF(AIYY.LE.DLIMIT) ILIN=1
IF(AIZZ.LE.DLIMIT) ILIN=1
WRITE(6,4)
CALL PRINT(T,6,3,6)
NVIB=N3N-6+ILIN
C
C THE CALCULATION OF PRINCIPAL MOMENTS OF INERTIA
CALL RSP(3,3,6,T,PT,1,EG,FV1,FV2)
CALL ZERO(ROTAA,3)
CALL ZERO(ROTGC,3)
CALL ZERO(ROTCM,3)
CALL ZERO(ROTMH,3)
WRITE(6,5)
II=0
DO 106 I=1,3
PA=A00
PB=A00
PC=A00
IF(DABS(PT(I)).LE.DLIMIT) GO TO 201
II=II+1
PA=PT(I)*PARA
PB=CONST/PT(I)
PC=CYCL/PT(I)
ROTAA(II)=PT(I)
ROTGC(II)=PA
ROTCM(II)=PB
ROTMH(II)=PC
201 CONTINUE
WRITE(6,6) RAXIS(I),PT(I),PA,PB,PC
106 CONTINUE
WRITE(6,7)
CALL MATOUT(EG,3,3,3,3,6)
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C% ASSIGN THE ROTATIONAL TYPE OF THE ISOTOPIC MOLECULE %
C% ITOP IS A PARAMETER FOR ROTATIONAL TOP %
C% ITOP = 1 : DIATOMIC MOLECULE %
C% ITOP = 2 : LINEAR POLYATOMIC MOLECULE %
C% ITOP = 3 : ASYMMETRIC TOP %
C% ITOP = 4 : SYMMETRIC TOP (PROLATE) %
C% ITOP = 5 : SYMMETRIC TOP (OBLATE) %
C% ITOP = 6 : SPHERICAL TOP %
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ITOP=3
IF(DABS(ROTAA(2)-ROTAA(3)).LE.RLIMIT) ITOP=4
IF(DABS(ROTAA(1)-ROTAA(2)).LE.RLIMIT) ITOP=5
IF(DABS(ROTAA(1)-ROTAA(2)).LE.RLIMIT.AND.
* DABS(ROTAA(2)-ROTAA(3)).LE.RLIMIT) ITOP=6
IF(ILIN.NE.0) ITOP=2
IF(NATOM.EQ.2) ITOP=1
WRITE(6,8) RTYPE(ITOP)
C
C ASSIGN PRINCIPAL AXIS
IAXIS(1)=1
IAXIS(2)=2
IAXIS(3)=3
C
C FORM UX MATRIX
CALL ZERO(UX,N3N*N3N)
C
C%%%%%%%%%%%%%%%%%%
C% FOR DIATOMIC %
C%%%%%%%%%%%%%%%%%%
IF(ITOP.EQ.1) THEN
CALL ZERO(EG,9)
EG(1,1)=ONE
EG(2,2)=ONE
EG(3,3)=ONE
END IF
C
C FOR LINEAR POLYATOMIC
IF(ITOP.EQ.2) THEN
DO 108 I=1,3
TG=EG(I,1)
EG(I,1)=EG(I,3)
EG(I,3)=TG
108 CONTINUE
END IF
C
C%%%%%%%%%%%%%%%%%%%%%%%%
C% FOR ASYMMETRIC TOP %
C%%%%%%%%%%%%%%%%%%%%%%%%
C FOR PROLATE TOP RAY'S KAPPA = -1
C FOR OBLATE TOP RAY'S KAPPA = +1
IF(ITOP.EQ.3) THEN
BP=(ROTCM(3)-ROTCM(2))/(ROTCM(1)*TWO-ROTCM(2)-ROTCM(3))
BO=(ROTCM(1)-ROTCM(2))/(ROTCM(3)*TWO-ROTCM(2)-ROTCM(1))
RAYK=(ROTCM(2)*TWO-ROTCM(1)-ROTCM(3))/(ROTCM(1)-ROTCM(3))
WRITE(6,9) BP,BO,RAYK
C FOR AN ASYMMETRIC TOP THE COORDINATE SYSTEM WILL BE ROTATED
C TO SATISFY THE FOLLOWING RELATIONSHIPS
C I(X)<I(Y)<I(Z)
C I(A)<I(B)<I(C)
C A > B > C
C
IAXIS(1)=1
IAXIS(2)=2
IAXIS(3)=3
END IF
C
C%%%%%%%%%%%%%%%%%%%%%%%
C% FOR SYMMETRIC TOP %
C%%%%%%%%%%%%%%%%%%%%%%%
C FOR A SYMMETRIC TOP THE MOLECULAR AXIS IS ALONG THE Z AXIS
IF(ITOP.EQ.4.OR.ITOP.EQ.5) THEN
CALL ZERO(EG,9)
EG(1,1)=ONE
EG(2,2)=ONE
EG(3,3)=ONE
END IF
C FOR A PROLATE TOP
C I(X)<I(Y)=I(Z)
C I(A)<I(B)=I(C)
C A > B = C
IF(ITOP.EQ.4) THEN
IAXIS(1)=3
IAXIS(2)=2
IAXIS(3)=1
END IF
C
C FOR AN OBLATE TOP
C I(X)=I(Y)<I(Z)
C I(A)=I(B)<I(C)
C A = B > C
IF(ITOP.EQ.5) THEN
IAXIS(1)=1
IAXIS(2)=2
IAXIS(3)=3
END IF
C
C%%%%%%%%%%%%%%%%%%%%%%%
C% FOR SPHERICAL TOP %
C%%%%%%%%%%%%%%%%%%%%%%%
IF(ITOP.EQ.6) THEN
CALL ZERO(EG,9)
EG(1,1)=ONE
EG(2,2)=ONE
EG(3,3)=ONE
IAXIS(1)=1
IAXIS(2)=2
IAXIS(3)=3
END IF
C
C FORM THE UX MATRIX
DO 111 IABC=1,NATOM
IPOS=(IABC-1)*3
DO 110 I=1,3
DO 110 J=1,3
II=IPOS+I
JJ=IPOS+J
UX(II,JJ)=EG(I,J)
110 CONTINUE
111 CONTINUE
WRITE(6,10)
CALL MATOUT(UX,N3N,N3N,N3N,N3N,6)
C
C CHECK UNITARITY OF UX MATRIX
CALL MTXMPY(UX,N3N,UX,N3N,EE,N3N,EE,N3N,N3N,3)
WRITE(6,11)
CALL MATOUT(EE,N3N,N3N,N3N,N3N,6)
C
C THE CARTESIAN COORDINATES W.R.T. PRINCIPAL AXIS
WRITE(6,12)
DO 112 I=1,NATOM
XT=EG(1,1)*X(I)+EG(2,1)*Y(I)+EG(3,1)*Z(I)
YT=EG(1,2)*X(I)+EG(2,2)*Y(I)+EG(3,2)*Z(I)
ZT=EG(1,3)*X(I)+EG(2,3)*Y(I)+EG(3,3)*Z(I)
X(I)=XT
Y(I)=YT
Z(I)=ZT
XT=XT/A0
YT=YT/A0
ZT=ZT/A0
WRITE(6,3) I,XT,YT,ZT,W(I)
112 CONTINUE
WRITE(6,13)
DO 113 I=1,NATOM
113 WRITE(6,3) I,X(I),Y(I),Z(I),W(I)
C
RETURN
END
SUBROUTINE VIBFX(FX,ELX,ELXM,UX,EPX,EPXM,EE,VALU,FV1,FV2,FXM,DD,
1 CC)
C THE NORMAL COORDINATE ANALYSIS FOR CARTESIAN COORDINATE SYSTEM
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION FX(N3N,N3N),ELX(N3N,N3N),ELXM(N3N,N3N)
DIMENSION UX(N3N,N3N),EPX(N3N,N3N),EPXM(N3N,N3N)
DIMENSION EE(N3N,N3N),DD(N3N,N3N),CC(1)
DIMENSION VALU(N3N),FV1(N3N),FV2(N3N),FXM(NATRI)
DIMENSION IATOMX(10)
COMMON/VIB101/NATOM,N3N,NATRI,ILIN,NVIB
COMMON/VIB107/ITOP,IVSGL,IVDBL,IVTRP
COMMON/VIB108/IAXIS(3),NDEG(150),NDAB(150,5),IMAG(150)
COMMON/VIB201/CHG(50),X(50),Y(50),Z(50),W(50)
COMMON/VIB203/IOFF(150),IPRNT
COMMON/VIB204/SQM(150),ROOT(150),FREQ(150)
DATA A00,ONE / 0.0D+00 , 1.0D+00 /
DATA DLIMIT / 1.0D-01 /
DATA ELIMIT / 1.0D-04 /
DATA WLIMIT / 3.0D+00 /
DATA XLIMIT / 1.0D-05 /
1 FORMAT(//,2X,' THE NORMAL COORDINATE ANALYSIS BY FX MATRIX METHOD'
1/)
2 FORMAT(//,2X,' SECULAR EQUATIONS'/2X,' FXM MATRIX'/)
3 FORMAT(//,2X,' *****************************'/
1 2X,' ***VIBRATIONAL FREQUENCIES***'/
2 2X,' *****************************')
4 FORMAT(//,2X,' LX MATRIX'/)
5 FORMAT(//,2X,' LX * LX(T) = M(-1)'/)
6 FORMAT(//,2X,' LAMDA = LX(T) * FX * LX'/)
7 FORMAT(//,2X,' LX'' MATRIX'/)
8 FORMAT(//,2X,' LX'' * LX''(T) = M(-1)'/)
9 FORMAT(//,2X,' LXM'' MATRIX'/)
10 FORMAT(//,2X,' LXM'' * LXM''(T) = E'/)
11 FORMAT(//,2X,' ANALYSIS OF VIBRATIONAL MODES'//
1 2X,' NO.OF NON-DEGENERATE MODES = ',I5/
1 2X,' NO.OF DOUBLY-DEGENERATE MODES = ',I5/
2 2X,' NO.OF TRIPLY-DEGENERATE MODES = ',I5/
3 2X,' NO.OF TOTAL VIBRATIONAL MODES = ',I5/)
12 FORMAT(//,2X,' ROTATED LXM'' MATRIX'/)
C
IF(IPRNT.NE.0)
*WRITE(6,1)
C
C THE FORMATION OF THE SECULAR EQUATION
DO 101 I=1,NATOM
101 SQM(I)=DSQRT(W(I))
DO 102 I=1,N3N
II=I/3
IF(I.GT.3*II) II=II+1
DO 102 J=I,N3N
JJ=J/3
IF(J.GT.3*JJ) JJ=JJ+1
IJ=I+IOFF(J)
FXM(IJ)=(ONE/SQM(II))*FX(I,J)*(ONE/SQM(JJ))
102 CONTINUE
WRITE(6,2)
CALL PRINT(FXM,NATRI,N3N,6)
C
C THE NORMAL COORDINATE ANALYSIS
WRITE(6,3)
IF(ITOP.EQ.2) THEN
NAA=IOFF(NATOM+1)
IC1=1
IC2=IC1+NAA
IC3=IC2+NATOM
IC4=IC3+NATOM*NATOM
ICMAX=IC4+N3N
C
C................... FXMX VALUX ELXMX...
CALL NORMLI(ELXM,EE,VALU,FV1,FV2,FXM,CC(IC1),CC(IC2),CC(IC3),
C...................NORD........
1 CC(IC4),NAA)
ELSE
CALL NORMFX(ELXM,EE,VALU,FV1,FV2,FXM)
END IF
C
C THE CALCULATION OF LX MATRIX
DO 103 I=1,N3N
II=I/3
IF(I.GT.3*II) II=II+1
DO 103 J=1,N3N
ELX(I,J)=(ONE/SQM(II))*ELXM(I,J)
103 CONTINUE
WRITE(6,4)
CALL FRQOUT(ELX,FREQ,N3N,N3N,N3N,N3N,6)
IF(IPRNT.LE.2) GO TO 201