Revision 5 src/Read_geom.f90
Read_geom.f90 (revision 5)  

1  1 
SUBROUTINE Read_Geom(input) 
2  2  
3 
Use Path_module 

3 
Use VarTypes 

4 
Use Path_module, only : NGeomI 

4  5 
Use Io_module 
5  6  
6  7 
IMPLICIT NONE 
...  ...  
42  43  
43  44 
SELECT CASE(Input) 
44  45 
CASE ('XYZ','CART') 
45 
DO I=1,NGeomI 

46 
IF (DEBUG) WRITE(*,*) "Reading Geom :",I 

47 
READ(IOIN,*) NAtp 

48 
if (NAtp.NE.Nat) THEN 

49 
IF (I==1) THEN 

50 
WRITE(IOOUT,*) 'WARNING Number of atoms not consistent between NAMELIST &path and First geom' 

51 
WRITE(IOOUT,*) "Using:",Natp 

52 
DEALLOCATE(XyzGeomI, AtName) 

53 
Nat=Natp 

54 
ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt)) 

55 
ELSE 

56 
WRITE(IOOUT,*) 'Number of atoms not consistent between geometries. STOP' 

57 
STOP 

58 
END IF 

59 
END IF 

60 
READ(IOIN,'(A)') Line 

61 
DO J=1,NAt 

62 
READ(IOIN,*) AtName(J),XyzGeomI(I,1:3,J) 

63 
END DO 

64 
If (Debug) THEN 

65 
WRITE(*,*) "Geom ",I 

66 
DO J=1,NAt 

67 
WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(J),XyzGeomI(I,1:3,J) 

68 
END DO 

69 
END IF 

70 
END DO 

46 
Call ReadGeom_cart 

71  47 
CASE ('TURBOMOLE') 
72 
DO I=1,NGeomI 

73 
IF (DEBUG) WRITE(*,*) "Reading Geom :",I 

74 
JStart=1 

75 
! read the $coord line, or $userdefined, or $end 

76 
Line='$ok' 

77 
DO WHILE (Line(1:1).EQ."$") 

78 
READ(IOIN,'(A)') Line 

79 
Line=AdjustL(Line) 

80 
END DO 

81 
J=1 

82 
READ(Line,*) XyzGeomI(I,1:3,J),AtName(J) 

83 
DO J=2,NAt 

84 
READ(IOIN,*) XyzGeomI(I,1:3,J),AtName(J) 

85 
END DO 

86 
XyzGeomI(I,:,:)= XyzGeomI(I,:,:)*a0 

87 
If (Debug) THEN 

88 
WRITE(*,*) "Geom ",I 

89 
DO J=1,NAt 

90 
WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(J),XyzGeomI(I,1:3,J) 

91 
END DO 

92 
END IF 

93 
END DO 

48 
Call ReadGeom_turbomole 

94  49 
CASE ('VASP') 
95 
ALLOCATE(FFF(3,nat),FFFl(3,Nat),FFFt(3,Nat)) 

96 
! First geometry is a bit special for VASP as we have to set 

97 
! many things 

98 
IGeom=1 

99 
IF (DEBUG) WRITE(*,*) "Reading Geom :",IGeom 

100 
READ(IOIN,'(A)') Vasp_Title 

101 
READ(IOIN,*) Vasp_param 

102  
103 
READ(IOIN,*) Lat_a 

104 
READ(IOIN,*) Lat_b 

105 
READ(IOIN,*) Lat_c 

106  
107 
Lat_a=Lat_a*Vasp_param 

108 
Lat_b=Lat_b*Vasp_param 

109 
Lat_c=Lat_c*Vasp_param 

110  
111 
ALLOCATE(NbAtType(nat)) 

112 
READ(IOIN,'(A)') Vasp_types 

113 
! First, how many different types ? 

114 
NbAtType=0 

115 
READ(Vasp_types,*,END=999) NbAtType 

116 
999 CONTINUE 

117 
NbType=0 

118 
NatP=0 

119 
DO WHILE (NbAtType(NbType+1).NE.0) 

120 
NbType=NbType+1 

121 
Natp=Natp+NbAtType(NbType) 

122 
END DO 

123  
124  
125 
if (NAtp.NE.Nat) THEN 

126 
WRITE(IOOUT,*) 'WARNING Number of atoms not consistent between NAMELIST &path and First geom' 

127 
WRITE(IOOUT,*) "Using:",Natp 

128 
DEALLOCATE(XyzGeomI, AtName) 

129 
Nat=Natp 

130 
ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt)) 

131 
END if 

132  
133 
! Do we know the atom types ? 

134 
IF (AtTypes(1).EQ.' ') THEN 

135 
! user has not provided atom types... we have to find them by ourselves 

136 
! by looking into the POTCAR file... 

137 
INQUIRE(File="POTCAR", EXIST=Tchk) 

138 
IF (.NOT.Tchk) THEN 

139 
WRITE(*,*) "ERROR! No AtTypes provided, and POTCAR file not found" 

140 
STOP 

141 
END IF 

142 
OPEN(IOTMP,File="POTCAR") 

143 
DO I=1,NbType 

144 
Line='Empty' 

145 
DO WHILE (Line(1:5).NE.'TITEL') 

146 
READ(IOTMP,'(A)') Line 

147 
Line=AdjustL(Line) 

148 
END DO 

149 
Line=adjustl(Line(9:)) 

150 
Idx=index(Line," ") 

151 
Line=adjustl(Line(Idx+1:)) 

152 
AtTypes(I)=Line(1:2) 

153 
END DO 

154 
if (debug) WRITE(*,'(A,100(1X,A2))') "ReadG:VASP AtTypes",AtTypes(1:NbType) 

155 
CLOSE(IOTMP) 

156  
157 
ELSE !AtTypes(1).EQ.' ' 

158 
! user has provided atom types 

159 
NbTypeUser=0 

160 
DO WHILE (AtTypes(NbTypeUser+1).NE.' ') 

161 
NbTypeUser=NbTypeUser+1 

162 
END DO 

163 
IF (NbType.NE.NbTypeUser) THEN 

164 
WRITE(*,*) "ERROR Read_Geom : NbType in POSCAR do not match AtTypes" 

165 
STOP 

166 
END IF 

167 
END IF 

168  
169 
IAt=1 

170 
DO I=1,NbType 

171 
DO J=1,NbAtType(I) 

172 
AtName(Iat)=AtTypes(I) 

173 
Iat=Iat+1 

174 
END DO 

175 
END DO 

176 
DEALLOCATE(NbAtType) 

177  
178 
NbTypes=NbType 

179  
180 
READ(IOIN,'(A)') Vasp_comment 

181 
READ(IOIN,'(A)') Vasp_direct 

182  
183 
V_direct=adjustl(Vasp_direct) 

184 
Call upcase(V_direct) 

185  
186 
if (debug) WRITE(*,*) "V_direct=",Trim(V_Direct) 

187  
188 
DO I=1,Nat 

189 
READ(IOIN,*) XyzGeomI(IGeom,1:3,I),FFF(1:3,I) 

190 
DO J=1,3 

191 
FFF(J,I)=AdjustL(FFF(J,I)) 

192 
FFFl(j,i)=(FFF(j,i)(1:1).EQ.'T').OR.(FFF(j,i)(1:1).EQ.'t') 

193 
END DO 

194  
195 
END DO 

196 
IF (V_direct(1:6).EQ.'DIRECT') THEN 

197 
DO I=1,Nat 

198 
Xt=XyzGeomI(IGeom,1,I) 

199 
Yt=XyzGeomI(IGeom,2,I) 

200 
Zt=XyzGeomI(IGeom,3,I) 

201 
XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1) 

202 
XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2) 

203 
XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3) 

204 
END DO 

205  
206 
END IF 

207  
208 
ALLOCATE(X0_vasp(Nat),Y0_vasp(Nat),Z0_vasp(Nat)) 

209 
X0_vasp=XyzGeomI(IGeom,1,:) 

210 
Y0_vasp=XyzGeomI(IGeom,2,:) 

211 
Z0_vasp=XyzGeomI(IGeom,3,:) 

212  
213 
if (debug) THEN 

214 
WRITE(*,*) Nat 

215 
WRITE(*,*) "Geom ",IGeom,"/",NGeomI 

216 
DO I=1,Nat 

217 
WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(I),XyzGeomI(IGeom,1:3,I) 

218 
END DO 

219 
END IF 

220  
221  
222 
DO IGeom=2,NGeomI 

223 
!!! PFL 16th March 2008 > 

224 
! Users might want to use different conditions for different geometries 

225 
! i.e. to mix Cartesian and Direct :p 

226 
! thus we check for each geom which condition is used. 

227 
! 

228  
229 
! Title 

230 
READ(IOIN,*) LineTmp 

231 
! if (debug) WRITE(*,*) "1 LineTmp:",Trim(LineTmp) 

232 
! Vasp_Param 

233 
READ(IOIN,*) VaspPar_loc 

234 
! if (debug) WRITE(*,*) "2 LineTmp:",Trim(LineTmp) 

235 
! Lattice parameters 

236 
READ(IOIN,*) Lata_loc 

237 
! if (debug) WRITE(*,*) "3 LineTmp:",Trim(LineTmp) 

238 
READ(IOIN,*) latb_loc 

239 
! if (debug) WRITE(*,*) "4 LineTmp:",Trim(LineTmp) 

240 
READ(IOIN,*) Latc_loc 

241 
! if (debug) WRITE(*,*) "5 LineTmp:",Trim(LineTmp) 

242 
! Attypes 

243 
READ(IOIN,*) LineTmp 

244 
! if (debug) WRITE(*,*) "6 LineTmp:",Trim(LineTmp) 

245 
! Comment 

246 
READ(IOIN,*) LineTmp,line 

247 
! if (debug) WRITE(*,*) "7 LineTmp:",Trim(LineTmp),Trim(Line) 

248 
! V_direct 

249 
READ(IOIN,*) Vdir_loc 

250 
Vdir_loc=adjustl(Vdir_loc) 

251 
Call upcase(Vdir_loc) 

252 
! if (debug) WRITE(*,*) "8 LineTmp:",Trim(LineTmp) 

253  
254 
Lata_loc=Lata_loc*VaspPar_locLat_a 

255 
if (dot_product(lata_loc,Lata_loc).GE.1e6) THEN 

256 
WRITE(*,*) "Lattice parameter a is different between geometries 1 and ",IGeom 

257 
WRITE(IOOUT,*) "Lattice parameter a is different between geometries 1 and ",IGeom 

258 
STOP 

259 
END IF 

260  
261 
Latb_loc=Latb_loc*VaspPar_locLat_b 

262 
if (dot_product(latb_loc,Latb_loc).GE.1e6) THEN 

263 
WRITE(*,*) "Lattice parameter b is different between geometries 1 and ",IGeom 

264 
WRITE(IOOUT,*) "Lattice parameter b is different between geometries 1 and ",IGeom 

265 
STOP 

266 
END IF 

267  
268 
Latc_loc=Latc_loc*VaspPar_locLat_c 

269 
if (dot_product(latc_loc,Latc_loc).GE.1e6) THEN 

270 
WRITE(*,*) "Lattice parameter c is different between geometries 1 and ",IGeom 

271 
WRITE(IOOUT,*) "Lattice parameter c is different between geometries 1 and ",IGeom 

272 
STOP 

273 
END IF 

274  
275 
DO I=1,Nat 

276 
READ(IOIN,*) XyzGeomI(IGeom,1:3,I),FFFt(1:3,I) 

277 
DO J=1,3 

278 
! A frozen atom in any geom is considered frozen in all geom 

279 
FFFt(j,I)=AdjustL(FFFt(j,i)) 

280 
FFFl(j,i)=FFFl(j,i).AND.(FFFt(j,i).EQ.'T') 

281 
END DO 

282 
END DO 

283  
284 
IF (Vdir_loc(1:6).EQ.'DIRECT') THEN 

285 
DO I=1,Nat 

286 
Xt=XyzGeomI(IGeom,1,I) 

287 
Yt=XyzGeomI(IGeom,2,I) 

288 
Zt=XyzGeomI(IGeom,3,I) 

289 
XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1) 

290 
XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2) 

291 
XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3) 

292 
END DO 

293 
END IF 

294  
295 
if (debug) THEN 

296 
WRITE(*,*) Nat 

297 
WRITE(*,*) "Geom ",IGeom,"/",NGeomI 

298 
DO I=1,Nat 

299 
WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(I),XyzGeomI(IGeom,1:3,I) 

300 
END DO 

301 
END IF 

302  
303 
END DO 

304  
305 
! PFL 23 Nov 2008 > 

306 
! This is done in Path.f90 

307 
! ! We put FFrozen in place 

308 
! Idx=1 

309 
! DO I=1,Nat 

310 
! TChk=.FALSE. 

311 
! DO J=1,3 

312 
! TChk=Tchk.OR.(.NOT.FFFl(j,i)) 

313 
! END DO 

314 
! IF (TChk) THEN 

315 
! Frozen(Idx)=I 

316 
! Idx=Idx+1 

317 
! END IF 

318 
! END DO 

319 
! Frozen(Idx)=0 

320 
! We put FFF in place because FFFl is a logical flag 

321 
DO I=1,Nat 

322 
DO J=1,3 

323 
if (FFFl(j,i)) FFF(j,i)="T" 

324 
END DO 

325 
END DO 

326  
327 
!< PFL 23 Nov 2008 

328  
329 
Deallocate(FFFl,FFFt) 

330  
50 
Call ReadGeom_vasp 

51 
CASE ('SIESTA') 

52 
Call ReadGeom_siesta 

331  53 
CASE Default 
332  54 
WRITe(*,*) 'Input=',trim(Input),' UNKNOWN. Stop' 
333  55 
STOP 
Also available in: Unified diff