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
|
SUBROUTINE GRID
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'SIZES'
************************************************************************
*
* GRID CALCULATES THE ENERGY-SURFACE RESULTING FROM VARIATION OF
* TWO COORDINATES. THE STEP-SIZE IS STEP1 AND STEP2, AND A 11
* BY 11 GRID OF POINTS IS GENERATED
*
* This subroutine is extensively modified by Manyin Yi, Aug 1989.
*
* New features:
* 1. The input geometry definition should always be the upper-left
* corner(smallest coordinates) instead of the middle point;
* 2. The starting point for calculation can be one of the four
* corners by setting "+/-" STEP1/2;
* 3. The grid size(max 23*23) is controlled by POINT1/2,
* if POINT1/2 is omitted, then a size of 11*11 is assumed;
* kwd MAX sets the max size;
* 4. The upper-left corner of the plotting grid always corresponds
* to the smallest coordinates, the lower-right corner to the
* largest, no matter where the calculation starts;
* 5. Restartable.
* 6. Write out UNIMAP irregular data UMP.DAT
*
************************************************************************
COMMON /GEOM / GEO(3,NUMATM), XCOORD(3,NUMATM)
COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
COMMON /GRADNT/ GRAD(MAXPAR),GNORM
COMMON /GRAVEC/ COSINE
COMMON /MESH / LATOM1, LPARA1, LATOM2, LPARA2
COMMON /GPARAM/ CURRT1,CURRT2
COMMON /IJLP / IJLP, ILP, JLP, JLP1, IONE
COMMON /SURF / SURF
COMMON /KEYWRD/ KEYWRD
COMMON /TITLES/ KOMENT, TITLE
CHARACTER KEYWRD*241, KOMENT*81, TITLE*81, GETNAM*80
LOGICAL RESTRT
DIMENSION GD(MAXPAR),XLAST(MAXPAR),MDFP(20),XDFP(20)
DIMENSION SURFAC(23,23)
DIMENSION SURF(23*23)
DIMENSION UMPX(23),UMPY(23),UMPZ(23*23)
C
STEP1=READA(KEYWRD,INDEX(KEYWRD,'STEP1')+6)
STEP2=READA(KEYWRD,INDEX(KEYWRD,'STEP2')+6)
NPTS1=11
NPTS2=11
IF (INDEX(KEYWRD,' MAX').NE.0) THEN
NPTS1=23
NPTS2=23
GOTO 10
ENDIF
IF (INDEX(KEYWRD,'POINT1').NE.0)
1 NPTS1=ABS(READA(KEYWRD,INDEX(KEYWRD,'POINT1')+7))
IF (INDEX(KEYWRD,'POINT2').NE.0)
1 NPTS2=ABS(READA(KEYWRD,INDEX(KEYWRD,'POINT2')+7))
10 RESTRT=(INDEX(KEYWRD,'RESTART').NE.0)
C
C THE TOP-LEFT VALUE OF THE FIRST AND SECOND DIMENSIONS ARE
C GEO(LPARA1,LATOM1) AND GEO(LPARA2,LATOM2)
C
UMPY(1)=GEO(LPARA1,LATOM1)
UMPX(1)=GEO(LPARA2,LATOM2)
DEGREE=180.D0/3.14159265359D0
IF(LPARA1.NE.1)STEP1=STEP1/DEGREE
IF(LPARA2.NE.1)STEP2=STEP2/DEGREE
C
C NOW SET THE STARTING POINT TO THE DESIRED CORNER
C
IF(STEP1.GT.0.0.AND.STEP2.GT.0.0) THEN
START1=GEO(LPARA1,LATOM1)
START2=GEO(LPARA2,LATOM2)
ENDIF
C BOTTOM-LEFT
IF(STEP1.LT.0.0.AND.STEP2.GT.0.0)THEN
START1=GEO(LPARA1,LATOM1)+(NPTS1-1)*(ABS(STEP1))
START2=GEO(LPARA2,LATOM2)
ENDIF
C TOP-RIGHT
IF(STEP1.GT.0.0.AND.STEP2.LT.0.0)THEN
START1=GEO(LPARA1,LATOM1)
START2=GEO(LPARA2,LATOM2)+ABS((NPTS2-1)*STEP2)
ENDIF
C BOTTOM-RIGHT
IF(STEP1.LT.0.0.AND.STEP2.LT.0.0)THEN
START1=GEO(LPARA1,LATOM1)+ABS((NPTS1-1)*STEP1)
START2=GEO(LPARA2,LATOM2)+ABS((NPTS2-1)*STEP2)
ENDIF
C
C NOW TO SWEEP THROUGH THE GRID OF POINTS LEFT TO RIGHT THEN RIGHT
C TO LEFT OR VISA VERSA. THIS SHOULD AVOID THE GEOMETRY OR SCF GETTING
C MESSED UP.
C
IF(LPARA1.NE.1) THEN
C1=DEGREE
ELSE
C1=1.D0
ENDIF
IF(LPARA2.NE.1) THEN
C2=DEGREE
ELSE
C2=1.D0
ENDIF
C THESE PARAMETERS NEED TO BE DUMPED IN '.RES'
CURRT1=START1
CURRT2=START2
IONE=-1
CPUTOT=0.0D0
IJLP=0
ILP=1
JLP=1
JLP1=1
SURF(1)=0.D0
C
IF (RESTRT) THEN
MDFP(9)=0
CALL DFPSAV(CPUTOT,XPARAM,GD,XLAST,ESCF,MDFP,XDFP)
ENDIF
C
GEO(LPARA1,LATOM1)=CURRT1
GEO(LPARA2,LATOM2)=CURRT2
DO 30 ILOOP=ILP,NPTS1
IONE=-IONE
DO 20 JLOOP=JLP,NPTS2
JLOOP1=0
IF(IONE.LT.0)JLOOP1=NPTS2+1
IF(RESTRT) THEN
JLOOP1=JLP1
IONE=-IONE
RESTRT=.FALSE.
ELSE
JLOOP1=JLOOP1+IONE
JLP1=JLOOP1
ENDIF
CPU1=SECOND()
CURRT1=GEO(LPARA1,LATOM1)
CURRT2=GEO(LPARA2,LATOM2)
CALL FLEPO(XPARAM, NVAR, ESCF)
CPU2=SECOND()
CPU3=CPU2-CPU1
CPUTOT=CPUTOT+CPU3
JLP=JLP+1
IJLP=IJLP+1
SURF(IJLP)=ESCF
WRITE(6,'(/'' FIRST VARIABLE SECOND VARIABLE
1 FUNCTION'')')
WRITE(6,'('' :'',F16.5,F16.5,F16.6)')GEO(LPARA1,LATOM1)*C1,
1 GEO(LPARA2,LATOM2)*C2,ESCF
CALL GEOUT(6)
GEO(LPARA2,LATOM2)=GEO(LPARA2,LATOM2)+STEP2*IONE
20 CONTINUE
GEO(LPARA1,LATOM1)=GEO(LPARA1,LATOM1)+STEP1
GEO(LPARA2,LATOM2)=GEO(LPARA2,LATOM2)-STEP2*IONE
ILP=ILP+1
JLP=1
30 CONTINUE
WRITE(6,'(/10X,''HORIZONTAL: VARYING SECOND PARAMETER,'',
1 /10X,''VERTICAL: VARYING FIRST PARAMETER'')')
WRITE(6,'(/10X,''WHOLE OF GRID, SUITABLE FOR PLOTTING'',//)')
C
C ARCHIVE
OPEN(UNIT=12,FILE=GETNAM('FOR012'),STATUS='UNKNOWN')
OPEN(UNIT=20,FILE=GETNAM('FOR020'),STATUS='NEW',ERR=31)
GOTO 32
31 OPEN(UNIT=20,FILE=GETNAM('FOR020'),STATUS='OLD')
32 CONTINUE
WRITE(12,40)
CALL WRTTXT(12)
40 FORMAT(' ARCHIVE FILE FOR GRID CALCULATION'/'GRID OF HEATS'/)
WRITE(12,'(/'' TOTAL CPU TIME IN FLEPO : '',F10.3/)') CPUTOT
C
C WRITE OUT THE GRIDS
IONE=1
ILOOP=1
JLOOP1=1
DO 50 IJ=1,NPTS1*NPTS2
SURFAC(JLOOP1,ILOOP)=SURF(IJ)
N=IJ-(IJ/NPTS2)*NPTS2
IF (N.EQ.0) THEN
ILOOP=ILOOP+1
JLOOP1=JLOOP1+IONE
IONE=-IONE
ENDIF
JLOOP1=JLOOP1+IONE
50 CONTINUE
C
DO 60 I=2,NPTS1
60 UMPY(I)=UMPY(1)+(I-1)*ABS(STEP1)
DO 70 I=2,NPTS2
70 UMPX(I)=UMPX(1)+(I-1)*ABS(STEP2)
N=0
IF(STEP1.GT.0.0.AND.STEP2.GT.0.0) THEN
DO 90 I=1,NPTS1
DO 80 J=1,NPTS2
N=N+1
80 UMPZ(N)=SURFAC(J,I)
WRITE(6,'(11F7.2)')(SURFAC(J,I),J=1,NPTS2)
90 WRITE(12,'(11F7.2)')(SURFAC(J,I),J=1,NPTS2)
ENDIF
IF(STEP1.LT.0.0.AND.STEP2.GT.0.0) THEN
DO 110 I=NPTS1,1,-1
DO 100 J=1,NPTS2
N=N+1
100 UMPZ(N)=SURFAC(J,I)
WRITE(6,'(11F7.2)')(SURFAC(J,I),J=1,NPTS2)
110 WRITE(12,'(11F7.2)')(SURFAC(J,I),J=1,NPTS2)
ENDIF
IF(STEP1.GT.0.0.AND.STEP2.LT.0.0) THEN
DO 130 I=1,NPTS1
DO 120 J=NPTS2,1,-1
N=N+1
120 UMPZ(N)=SURFAC(J,I)
WRITE(6,'(11F7.2)')(SURFAC(J,I),J=NPTS2,1,-1)
130 WRITE(12,'(11F7.2)')(SURFAC(J,I),J=NPTS2,1,-1)
ENDIF
IF(STEP1.LT.0.0.AND.STEP2.LT.0.0) THEN
DO 150 I=NPTS1,1,-1
DO 140 J=NPTS2,1,-1
N=N+1
140 UMPZ(N)=SURFAC(J,I)
WRITE(6,'(11F7.2)')(SURFAC(J,I),J=NPTS2,1,-1)
150 WRITE(12,'(11F7.2)')(SURFAC(J,I),J=NPTS2,1,-1)
ENDIF
DO 160 I=0,NPTS1-1
DO 160 J=1,NPTS2
N=I*NPTS2+J
160 WRITE(20,'(3(1X,F8.3))')UMPX(J),UMPY(I+1),UMPZ(N)
CLOSE(20)
END
|
