root / src / Read_geom.f90 @ 2
Historique | Voir | Annoter | Télécharger (9,8 ko)
1 | 1 | equemene | SUBROUTINE Read_Geom(input) |
---|---|---|---|
2 | 1 | equemene | |
3 | 1 | equemene | Use Path_module |
4 | 1 | equemene | Use Io_module |
5 | 1 | equemene | |
6 | 1 | equemene | IMPLICIT NONE |
7 | 1 | equemene | |
8 | 1 | equemene | |
9 | 1 | equemene | CHARACTER(32), INTENT(IN) :: input |
10 | 1 | equemene | CHARACTER(132) :: LineTmp,Line |
11 | 1 | equemene | |
12 | 1 | equemene | INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na |
13 | 1 | equemene | INTEGER(KINT) :: NbType, NbTypeUser |
14 | 1 | equemene | INTEGER(KINT) :: I, J, Iat,NAtP,Idx,JStart |
15 | 1 | equemene | INTEGER(KINT) :: IGeom |
16 | 1 | equemene | REAL(KREAL) :: Xt,Yt,Zt |
17 | 1 | equemene | LOGICAL :: Debug, TChk |
18 | 1 | equemene | CHARACTER(4), ALLOCATABLE :: FFFt(:,:) ! 3,Na |
19 | 1 | equemene | LOGICAL, ALLOCATABLE :: FFFl(:,:) ! 3,Na |
20 | 1 | equemene | |
21 | 1 | equemene | !!! for VASP inputs |
22 | 1 | equemene | ! latx_loc used to check that the lattice parameters are the same for all geometries |
23 | 1 | equemene | REAL(KREAL) :: lata_loc(3), latb_loc(3), latc_loc(3) |
24 | 1 | equemene | ! VaspPar_loc is the local value of the vasp parametre |
25 | 1 | equemene | REAL(KREAL) :: VaspPar_loc |
26 | 1 | equemene | ! Vdir_loc is used to convert the read geometry to Cartesian. |
27 | 1 | equemene | ! V_direct is set by the first geometry. |
28 | 1 | equemene | CHARACTER(LCHARS) :: Vdir_loc |
29 | 1 | equemene | |
30 | 1 | equemene | |
31 | 1 | equemene | INTERFACE |
32 | 1 | equemene | function valid(string) result (isValid) |
33 | 1 | equemene | CHARACTER(*), intent(in) :: string |
34 | 1 | equemene | logical :: isValid |
35 | 1 | equemene | END function VALID |
36 | 1 | equemene | END INTERFACE |
37 | 1 | equemene | |
38 | 1 | equemene | debug=valid('Read_geom') |
39 | 1 | equemene | if (debug) Call Header("Entering Read_Geom") |
40 | 1 | equemene | if (debug) WRITE(*,*) "Input:",Trim(Input) |
41 | 1 | equemene | |
42 | 1 | equemene | SELECT CASE(Input) |
43 | 1 | equemene | CASE ('XYZ','CART') |
44 | 1 | equemene | DO I=1,NGeomI |
45 | 1 | equemene | IF (DEBUG) WRITE(*,*) "Reading Geom :",I |
46 | 1 | equemene | READ(IOIN,*) NAtp |
47 | 1 | equemene | if (NAtp.NE.Nat) THEN |
48 | 1 | equemene | IF (I==1) THEN |
49 | 1 | equemene | WRITE(IOOUT,*) 'WARNING Number of atoms not consistent between NAMELIST &path and First geom' |
50 | 1 | equemene | WRITE(IOOUT,*) "Using:",Natp |
51 | 1 | equemene | DEALLOCATE(XyzGeomI, AtName) |
52 | 1 | equemene | Nat=Natp |
53 | 1 | equemene | ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt)) |
54 | 1 | equemene | ELSE |
55 | 1 | equemene | WRITE(IOOUT,*) 'Number of atoms not consistent between geometries. STOP' |
56 | 1 | equemene | STOP |
57 | 1 | equemene | END IF |
58 | 1 | equemene | END IF |
59 | 1 | equemene | READ(IOIN,'(A)') Line |
60 | 1 | equemene | DO J=1,NAt |
61 | 1 | equemene | READ(IOIN,*) AtName(J),XyzGeomI(I,1:3,J) |
62 | 1 | equemene | END DO |
63 | 1 | equemene | If (Debug) THEN |
64 | 1 | equemene | WRITE(*,*) "Geom ",I |
65 | 1 | equemene | DO J=1,NAt |
66 | 1 | equemene | WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(J),XyzGeomI(I,1:3,J) |
67 | 1 | equemene | END DO |
68 | 1 | equemene | END IF |
69 | 1 | equemene | END DO |
70 | 1 | equemene | CASE ('TURBOMOLE') |
71 | 1 | equemene | DO I=1,NGeomI |
72 | 1 | equemene | IF (DEBUG) WRITE(*,*) "Reading Geom :",I |
73 | 1 | equemene | JStart=1 |
74 | 1 | equemene | ! read the $coord line, or $user-defined, or $end |
75 | 1 | equemene | Line='$ok' |
76 | 1 | equemene | DO WHILE (Line(1:1).EQ."$") |
77 | 1 | equemene | READ(IOIN,'(A)') Line |
78 | 1 | equemene | Line=AdjustL(Line) |
79 | 1 | equemene | END DO |
80 | 1 | equemene | J=1 |
81 | 1 | equemene | READ(Line,*) XyzGeomI(I,1:3,J),AtName(J) |
82 | 1 | equemene | DO J=2,NAt |
83 | 1 | equemene | READ(IOIN,*) XyzGeomI(I,1:3,J),AtName(J) |
84 | 1 | equemene | END DO |
85 | 1 | equemene | XyzGeomI(I,:,:)= XyzGeomI(I,:,:)*a0 |
86 | 1 | equemene | If (Debug) THEN |
87 | 1 | equemene | WRITE(*,*) "Geom ",I |
88 | 1 | equemene | DO J=1,NAt |
89 | 1 | equemene | WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(J),XyzGeomI(I,1:3,J) |
90 | 1 | equemene | END DO |
91 | 1 | equemene | END IF |
92 | 1 | equemene | END DO |
93 | 1 | equemene | CASE ('VASP') |
94 | 1 | equemene | ALLOCATE(FFF(3,nat),FFFl(3,Nat),FFFt(3,Nat)) |
95 | 1 | equemene | ! First geometry is a bit special for VASP as we have to set |
96 | 1 | equemene | ! many things |
97 | 1 | equemene | IGeom=1 |
98 | 1 | equemene | IF (DEBUG) WRITE(*,*) "Reading Geom :",IGeom |
99 | 1 | equemene | READ(IOIN,'(A)') Vasp_Title |
100 | 1 | equemene | READ(IOIN,*) Vasp_param |
101 | 1 | equemene | |
102 | 1 | equemene | READ(IOIN,*) Lat_a |
103 | 1 | equemene | READ(IOIN,*) Lat_b |
104 | 1 | equemene | READ(IOIN,*) Lat_c |
105 | 1 | equemene | |
106 | 1 | equemene | Lat_a=Lat_a*Vasp_param |
107 | 1 | equemene | Lat_b=Lat_b*Vasp_param |
108 | 1 | equemene | Lat_c=Lat_c*Vasp_param |
109 | 1 | equemene | |
110 | 1 | equemene | ALLOCATE(NbAtType(nat)) |
111 | 1 | equemene | READ(IOIN,'(A)') Vasp_types |
112 | 1 | equemene | ! First, how many different types ? |
113 | 1 | equemene | NbAtType=0 |
114 | 1 | equemene | READ(Vasp_types,*,END=999) NbAtType |
115 | 1 | equemene | 999 CONTINUE |
116 | 1 | equemene | NbType=0 |
117 | 1 | equemene | NatP=0 |
118 | 1 | equemene | DO WHILE (NbAtType(NbType+1).NE.0) |
119 | 1 | equemene | NbType=NbType+1 |
120 | 1 | equemene | Natp=Natp+NbAtType(NbType) |
121 | 1 | equemene | END DO |
122 | 1 | equemene | |
123 | 1 | equemene | |
124 | 1 | equemene | if (NAtp.NE.Nat) THEN |
125 | 1 | equemene | WRITE(IOOUT,*) 'WARNING Number of atoms not consistent between NAMELIST &path and First geom' |
126 | 1 | equemene | WRITE(IOOUT,*) "Using:",Natp |
127 | 1 | equemene | DEALLOCATE(XyzGeomI, AtName) |
128 | 1 | equemene | Nat=Natp |
129 | 1 | equemene | ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt)) |
130 | 1 | equemene | END if |
131 | 1 | equemene | |
132 | 1 | equemene | ! Do we know the atom types ? |
133 | 1 | equemene | IF (AtTypes(1).EQ.' ') THEN |
134 | 1 | equemene | ! user has not provided atom types... we have to find them by ourselves |
135 | 1 | equemene | ! by looking into the POTCAR file... |
136 | 1 | equemene | INQUIRE(File="POTCAR", EXIST=Tchk) |
137 | 1 | equemene | IF (.NOT.Tchk) THEN |
138 | 1 | equemene | WRITE(*,*) "ERROR! No AtTypes provided, and POTCAR file not found" |
139 | 1 | equemene | STOP |
140 | 1 | equemene | END IF |
141 | 1 | equemene | OPEN(IOTMP,File="POTCAR") |
142 | 1 | equemene | DO I=1,NbType |
143 | 1 | equemene | Line='Empty' |
144 | 1 | equemene | DO WHILE (Line(1:5).NE.'TITEL') |
145 | 1 | equemene | READ(IOTMP,'(A)') Line |
146 | 1 | equemene | Line=AdjustL(Line) |
147 | 1 | equemene | END DO |
148 | 1 | equemene | Line=adjustl(Line(9:)) |
149 | 1 | equemene | Idx=index(Line," ") |
150 | 1 | equemene | Line=adjustl(Line(Idx+1:)) |
151 | 1 | equemene | AtTypes(I)=Line(1:2) |
152 | 1 | equemene | END DO |
153 | 1 | equemene | if (debug) WRITE(*,'(A,100(1X,A2))') "ReadG:VASP AtTypes",AtTypes(1:NbType) |
154 | 1 | equemene | CLOSE(IOTMP) |
155 | 1 | equemene | |
156 | 1 | equemene | ELSE !AtTypes(1).EQ.' ' |
157 | 1 | equemene | ! user has provided atom types |
158 | 1 | equemene | NbTypeUser=0 |
159 | 1 | equemene | DO WHILE (AtTypes(NbTypeUser+1).NE.' ') |
160 | 1 | equemene | NbTypeUser=NbTypeUser+1 |
161 | 1 | equemene | END DO |
162 | 1 | equemene | IF (NbType.NE.NbTypeUser) THEN |
163 | 1 | equemene | WRITE(*,*) "ERROR Read_Geom : NbType in POSCAR do not match AtTypes" |
164 | 1 | equemene | STOP |
165 | 1 | equemene | END IF |
166 | 1 | equemene | END IF |
167 | 1 | equemene | |
168 | 1 | equemene | IAt=1 |
169 | 1 | equemene | DO I=1,NbType |
170 | 1 | equemene | DO J=1,NbAtType(I) |
171 | 1 | equemene | AtName(Iat)=AtTypes(I) |
172 | 1 | equemene | Iat=Iat+1 |
173 | 1 | equemene | END DO |
174 | 1 | equemene | END DO |
175 | 1 | equemene | DEALLOCATE(NbAtType) |
176 | 1 | equemene | |
177 | 1 | equemene | NbTypes=NbType |
178 | 1 | equemene | |
179 | 1 | equemene | READ(IOIN,'(A)') Vasp_comment |
180 | 1 | equemene | READ(IOIN,'(A)') Vasp_direct |
181 | 1 | equemene | |
182 | 1 | equemene | V_direct=adjustl(Vasp_direct) |
183 | 1 | equemene | Call upcase(V_direct) |
184 | 1 | equemene | |
185 | 1 | equemene | if (debug) WRITE(*,*) "V_direct=",Trim(V_Direct) |
186 | 1 | equemene | |
187 | 1 | equemene | DO I=1,Nat |
188 | 1 | equemene | READ(IOIN,*) XyzGeomI(IGeom,1:3,I),FFF(1:3,I) |
189 | 1 | equemene | DO J=1,3 |
190 | 1 | equemene | FFF(J,I)=AdjustL(FFF(J,I)) |
191 | 1 | equemene | FFFl(j,i)=(FFF(j,i)(1:1).EQ.'T') |
192 | 1 | equemene | END DO |
193 | 1 | equemene | |
194 | 1 | equemene | END DO |
195 | 1 | equemene | IF (V_direct(1:6).EQ.'DIRECT') THEN |
196 | 1 | equemene | DO I=1,Nat |
197 | 1 | equemene | Xt=XyzGeomI(IGeom,1,I) |
198 | 1 | equemene | Yt=XyzGeomI(IGeom,2,I) |
199 | 1 | equemene | Zt=XyzGeomI(IGeom,3,I) |
200 | 1 | equemene | XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1) |
201 | 1 | equemene | XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2) |
202 | 1 | equemene | XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3) |
203 | 1 | equemene | END DO |
204 | 1 | equemene | |
205 | 1 | equemene | END IF |
206 | 1 | equemene | |
207 | 1 | equemene | ALLOCATE(X0_vasp(Nat),Y0_vasp(Nat),Z0_vasp(Nat)) |
208 | 1 | equemene | X0_vasp=XyzGeomI(IGeom,1,:) |
209 | 1 | equemene | Y0_vasp=XyzGeomI(IGeom,2,:) |
210 | 1 | equemene | Z0_vasp=XyzGeomI(IGeom,3,:) |
211 | 1 | equemene | |
212 | 1 | equemene | if (debug) THEN |
213 | 1 | equemene | WRITE(*,*) Nat |
214 | 1 | equemene | WRITE(*,*) "Geom ",IGeom,"/",NGeomI |
215 | 1 | equemene | DO I=1,Nat |
216 | 1 | equemene | WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(I),XyzGeomI(IGeom,1:3,I) |
217 | 1 | equemene | END DO |
218 | 1 | equemene | END IF |
219 | 1 | equemene | |
220 | 1 | equemene | |
221 | 1 | equemene | DO IGeom=2,NGeomI |
222 | 1 | equemene | !!! PFL 16th March 2008 -> |
223 | 1 | equemene | ! Users might want to use different conditions for different geometries |
224 | 1 | equemene | ! i.e. to mix Cartesian and Direct :-p |
225 | 1 | equemene | ! thus we check for each geom which condition is used. |
226 | 1 | equemene | ! |
227 | 1 | equemene | |
228 | 1 | equemene | ! Title |
229 | 1 | equemene | READ(IOIN,*) LineTmp |
230 | 1 | equemene | ! if (debug) WRITE(*,*) "1 LineTmp:",Trim(LineTmp) |
231 | 1 | equemene | ! Vasp_Param |
232 | 1 | equemene | READ(IOIN,*) VaspPar_loc |
233 | 1 | equemene | ! if (debug) WRITE(*,*) "2 LineTmp:",Trim(LineTmp) |
234 | 1 | equemene | ! Lattice parameters |
235 | 1 | equemene | READ(IOIN,*) Lata_loc |
236 | 1 | equemene | ! if (debug) WRITE(*,*) "3 LineTmp:",Trim(LineTmp) |
237 | 1 | equemene | READ(IOIN,*) latb_loc |
238 | 1 | equemene | ! if (debug) WRITE(*,*) "4 LineTmp:",Trim(LineTmp) |
239 | 1 | equemene | READ(IOIN,*) Latc_loc |
240 | 1 | equemene | ! if (debug) WRITE(*,*) "5 LineTmp:",Trim(LineTmp) |
241 | 1 | equemene | ! Attypes |
242 | 1 | equemene | READ(IOIN,*) LineTmp |
243 | 1 | equemene | ! if (debug) WRITE(*,*) "6 LineTmp:",Trim(LineTmp) |
244 | 1 | equemene | ! Comment |
245 | 1 | equemene | READ(IOIN,*) LineTmp,line |
246 | 1 | equemene | ! if (debug) WRITE(*,*) "7 LineTmp:",Trim(LineTmp),Trim(Line) |
247 | 1 | equemene | ! V_direct |
248 | 1 | equemene | READ(IOIN,*) Vdir_loc |
249 | 1 | equemene | Vdir_loc=adjustl(Vdir_loc) |
250 | 1 | equemene | Call upcase(Vdir_loc) |
251 | 1 | equemene | ! if (debug) WRITE(*,*) "8 LineTmp:",Trim(LineTmp) |
252 | 1 | equemene | |
253 | 1 | equemene | Lata_loc=Lata_loc*VaspPar_loc-Lat_a |
254 | 1 | equemene | if (dot_product(lata_loc,Lata_loc).GE.1e-6) THEN |
255 | 1 | equemene | WRITE(*,*) "Lattice parameter a is different between geometries 1 and ",IGeom |
256 | 1 | equemene | WRITE(IOOUT,*) "Lattice parameter a is different between geometries 1 and ",IGeom |
257 | 1 | equemene | STOP |
258 | 1 | equemene | END IF |
259 | 1 | equemene | |
260 | 1 | equemene | Latb_loc=Latb_loc*VaspPar_loc-Lat_b |
261 | 1 | equemene | if (dot_product(latb_loc,Latb_loc).GE.1e-6) THEN |
262 | 1 | equemene | WRITE(*,*) "Lattice parameter b is different between geometries 1 and ",IGeom |
263 | 1 | equemene | WRITE(IOOUT,*) "Lattice parameter b is different between geometries 1 and ",IGeom |
264 | 1 | equemene | STOP |
265 | 1 | equemene | END IF |
266 | 1 | equemene | |
267 | 1 | equemene | Latc_loc=Latc_loc*VaspPar_loc-Lat_c |
268 | 1 | equemene | if (dot_product(latc_loc,Latc_loc).GE.1e-6) THEN |
269 | 1 | equemene | WRITE(*,*) "Lattice parameter c is different between geometries 1 and ",IGeom |
270 | 1 | equemene | WRITE(IOOUT,*) "Lattice parameter c is different between geometries 1 and ",IGeom |
271 | 1 | equemene | STOP |
272 | 1 | equemene | END IF |
273 | 1 | equemene | |
274 | 1 | equemene | DO I=1,Nat |
275 | 1 | equemene | READ(IOIN,*) XyzGeomI(IGeom,1:3,I),FFFt(1:3,I) |
276 | 1 | equemene | DO J=1,3 |
277 | 1 | equemene | ! A frozen atom in any geom is considered frozen in all geom |
278 | 1 | equemene | FFFt(j,I)=AdjustL(FFFt(j,i)) |
279 | 1 | equemene | FFFl(j,i)=FFFl(j,i).AND.(FFFt(j,i).EQ.'T') |
280 | 1 | equemene | END DO |
281 | 1 | equemene | END DO |
282 | 1 | equemene | |
283 | 1 | equemene | IF (Vdir_loc(1:6).EQ.'DIRECT') THEN |
284 | 1 | equemene | DO I=1,Nat |
285 | 1 | equemene | Xt=XyzGeomI(IGeom,1,I) |
286 | 1 | equemene | Yt=XyzGeomI(IGeom,2,I) |
287 | 1 | equemene | Zt=XyzGeomI(IGeom,3,I) |
288 | 1 | equemene | XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1) |
289 | 1 | equemene | XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2) |
290 | 1 | equemene | XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3) |
291 | 1 | equemene | END DO |
292 | 1 | equemene | END IF |
293 | 1 | equemene | |
294 | 1 | equemene | if (debug) THEN |
295 | 1 | equemene | WRITE(*,*) Nat |
296 | 1 | equemene | WRITE(*,*) "Geom ",IGeom,"/",NGeomI |
297 | 1 | equemene | DO I=1,Nat |
298 | 1 | equemene | WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(I),XyzGeomI(IGeom,1:3,I) |
299 | 1 | equemene | END DO |
300 | 1 | equemene | END IF |
301 | 1 | equemene | |
302 | 1 | equemene | END DO |
303 | 1 | equemene | |
304 | 1 | equemene | ! PFL 23 Nov 2008 -> |
305 | 1 | equemene | ! This is done in Path.f90 |
306 | 1 | equemene | ! ! We put FFrozen in place |
307 | 1 | equemene | ! Idx=1 |
308 | 1 | equemene | ! DO I=1,Nat |
309 | 1 | equemene | ! TChk=.FALSE. |
310 | 1 | equemene | ! DO J=1,3 |
311 | 1 | equemene | ! TChk=Tchk.OR.(.NOT.FFFl(j,i)) |
312 | 1 | equemene | ! END DO |
313 | 1 | equemene | ! IF (TChk) THEN |
314 | 1 | equemene | ! Frozen(Idx)=I |
315 | 1 | equemene | ! Idx=Idx+1 |
316 | 1 | equemene | ! END IF |
317 | 1 | equemene | ! END DO |
318 | 1 | equemene | ! Frozen(Idx)=0 |
319 | 1 | equemene | ! We put FFF in place because FFFl is a logical flag |
320 | 1 | equemene | DO I=1,Nat |
321 | 1 | equemene | DO J=1,3 |
322 | 1 | equemene | if (FFFl(j,i)) FFF(j,i)="T" |
323 | 1 | equemene | END DO |
324 | 1 | equemene | END DO |
325 | 1 | equemene | |
326 | 1 | equemene | !<- PFL 23 Nov 2008 |
327 | 1 | equemene | |
328 | 1 | equemene | Deallocate(FFFl,FFFt) |
329 | 1 | equemene | |
330 | 1 | equemene | CASE Default |
331 | 1 | equemene | WRITe(*,*) 'Input=',trim(Input),' UNKNOWN. Stop' |
332 | 1 | equemene | STOP |
333 | 1 | equemene | |
334 | 1 | equemene | END SELECT |
335 | 1 | equemene | |
336 | 1 | equemene | |
337 | 1 | equemene | if (debug) Call Header("Exiting Read_Geom") |
338 | 1 | equemene | |
339 | 1 | equemene | END SUBROUTINE Read_Geom |