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
|
SUBROUTINE LOCAL(C,MDIM,NOCC,EIG)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'SIZES'
DIMENSION C(MDIM,MDIM), EIG(MAXORB)
COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
2 NCLOSE,NOPEN,NDUMY,FRACT
C**********************************************************************
C
C LOCALISATION SUBROUTINE
C ON INPUT
C C = EIGENVECTORS IN AN MDIM*MDIM MATRIX
C NOCC = NUMBER OF FILLED LEVELS
C NORBS = NUMBER OF ORBITALS
C NUMAT = NUMBER OF ATOMS
C NLAST = INTEGER ARRAY OF ATOM ORBITAL COUNTERS
C NFIRST = INTEGER ARRAY OF ATOM ORBITAL COUNTERS
C
C SUBROUTINE MAXIMIZES (PSI)**4
C REFERENCE_
C A NEW RAPID METHOD FOR ORBITAL LOCALISATION, P.G. PERKINS AND
C J.J.P. STEWART, J.C.S. FARADAY (II) 77, 000, (1981).
C
C MODIFIED AND CORRECTED TO AVOID SIGMA-PI ORBITAL MIXING BY
C JUAN CARLOS PANIAGUA, UNIVERSITY OF BARCELONA, MAY 1983.
C
C**********************************************************************
COMMON /SCRACH/ COLD(MAXORB,MAXORB),XDUMY(MAXPAR**2-MAXORB*MAXORB)
DIMENSION EIG1(MAXORB),PSI1(MAXORB),PSI2(MAXORB),
1 CII(MAXORB), REFEIG(MAXORB),IEL(20)
SAVE ELEMNT
CHARACTER*2 ELEMNT(99)
DATA ELEMNT/'H','HE',
1 'LI','BE','B','C','N','O','F','NE',
2 'NA','MG','AL','SI','P','S','CL','AR',
3 'K','CA','SC','TI','V','CR','MN','FE','CO','NI','CU',
4 'ZN','GA','GE','AS','SE','BR','KR',
5 'RB','SR','Y','ZR','NB','MO','TC','RU','RH','PD','AG',
6 'CD','IN','SN','SB','TE','I','XE',
7 'CS','BA','LA','CE','PR','ND','PM','SM','EU','GD','TB','DY',
8 'HO','ER','TM','YB','LU','HF','TA','W','RE','OS','IR','PT',
9 'AU','HG','TL','PB','BI','PO','AT','RN',
1 'FR','RA','AC','TH','PA','U','NP','PU','AM','CM','BK','CF','XX'/
NITER=100
EPS=1.0D-7
DO 10 I=1,NORBS
REFEIG(I)=EIG(I)
DO 10 J=1,NORBS
10 COLD(I,J)=C(I,J)
ITER=0
20 CONTINUE
SUM=0.D0
ITER=ITER+1
DO 80 I=1,NOCC
DO 70 J=1,NOCC
IF(J.EQ.I) GOTO 70
XIJJJ=0.0D0
XJIII=0.0D0
XIIII=0.0D0
XJJJJ=0.0D0
XIJIJ=0.0D0
XIIJJ=0.0D0
DO 30 K=1,NORBS
PSI1(K)=C(K,I)
30 PSI2(K)=C(K,J)
C NOW FOLLOWS THE RATE-DETERMINING STEP FOR THE CALCULATION
DO 50 K1=1,NUMAT
KL=NFIRST(K1)
KU=NLAST(K1)
DIJ=0.D0
DII=0.D0
DJJ=0.D0
C$DOIT ASIS
DO 40 K=KL,KU
DIJ=DIJ+PSI1(K)*PSI2(K)
DII=DII+PSI1(K)*PSI1(K)
DJJ=DJJ+PSI2(K)*PSI2(K)
40 CONTINUE
XIJJJ=XIJJJ+DIJ*DJJ
XJIII=XJIII+DIJ*DII
XIIII=XIIII+DII*DII
XJJJJ=XJJJJ+DJJ*DJJ
XIJIJ=XIJIJ+DIJ*DIJ
XIIJJ=XIIJJ+DII*DJJ
50 CONTINUE
AIJ=XIJIJ-(XIIII+XJJJJ-2.0D0*XIIJJ)/4.0D0
BIJ=XJIII-XIJJJ
CA=SQRT(AIJ*AIJ+BIJ*BIJ)
SA=AIJ+CA
IF(SA.LT.1.0D-14) GO TO 70
SUM=SUM+SA
CA=-AIJ/CA
CA=(1.0D0+SQRT((1.0D0+CA)/2.0D0))/2.0D0
IF((2.0D0*CA-1.0D0)*BIJ.LT.0.0D0)CA=1.0D0-CA
SA=SQRT(1.0D0-CA)
CA=SQRT(CA)
DO 60 K=1,NORBS
C(K,I)=CA*PSI1(K)+SA*PSI2(K)
60 C(K,J)=-SA*PSI1(K)+CA*PSI2(K)
70 CONTINUE
80 CONTINUE
SUM1=0.D0
DO 100 I=1,NOCC
DO 100 J=1,NUMAT
IL=NFIRST(J)
IU=NLAST(J)
X=0.0D0
C$DOIT ASIS
DO 90 K=IL,IU
90 X=X+C(K,I)**2
100 SUM1=SUM1+X*X
IF(SUM.GT.EPS.AND.ITER.LT.NITER) GO TO 20
WRITE(6,110)ITER,SUM1
110 FORMAT(/10X,'NUMBER OF ITERATIONS =',I4/
110X,'LOCALISATION VALUE =',F14.9,/)
WRITE(6,120)
120 FORMAT(3X,'NUMBER OF CENTERS',14X,'(COMPOSITION OF ORBITALS)'//)
DO 150 I=1,NOCC
SUM=0.D0
DO 140 J=1,NOCC
CO=0.D0
DO 130 K=1,NORBS
130 CO=CO+COLD(K,J)*C(K,I)
140 SUM=SUM+CO*CO*EIG(J)
150 EIG1(I)=SUM
DO 180 I=1,NOCC
X=100.D0
DO 160 J=I,NOCC
IF (X.LT.EIG1(J)) GOTO 160
X=EIG1(J)
I1=J
160 CONTINUE
EIG(I)=EIG1(I1)
X=EIG1(I1)
EIG1(I1)=EIG1(I)
EIG1(I)=X
DO 170 J=1,NORBS
X=C(J,I1)
C(J,I1)=C(J,I)
170 C(J,I)=X
180 CONTINUE
DO 250 I=1,NOCC
X=0.D0
DO 200 K1=1,NUMAT
KL=NFIRST(K1)
KU=NLAST(K1)
DII=0.D0
DO 190 K=KL,KU
190 DII=DII+C(K,I)**2
X=X+DII*DII
200 PSI1(K1)=DII*100.D0
X=1/X
DO 220 II=1,NUMAT
SUM=0.D0
DO 210 J=1,NUMAT
IF(PSI1(J).LT.SUM) GOTO 210
SUM=PSI1(J)
K=J
210 CONTINUE
PSI1(K)=0.D0
CII(II)=SUM
IEL(II)=K
IF(SUM.LT.1.D0) GOTO 230
220 CONTINUE
230 CONTINUE
II=II-1
WRITE(6,240)X,(ELEMNT(NAT(IEL(K))),IEL(K),CII(K),K=1,II)
240 FORMAT(F10.4,4(5(3X,A2,I3,F6.2),/10X))
250 CONTINUE
260 FORMAT(//20X,20H LOCALIZED ORBITALS ,//)
WRITE(6,260)
CALL MATOUT(C,EIG,NOCC,NORBS,MDIM)
270 FORMAT(10F12.6)
DO 280 I=1,NOCC
EIG(I)=REFEIG(I)
DO 280 J=1,NORBS
280 C(J,I)=COLD(J,I)
RETURN
END
|
