## root / src / ReadGeom_vasp.f90 @ 5

History | View | Annotate | Download (7.7 kB)

1 |
SUBROUTINE ReadGeom_vasp |
---|---|

2 | |

3 |
Use Path_module |

4 |
Use Io_module |

5 | |

6 |
IMPLICIT NONE |

7 | |

8 | |

9 |
CHARACTER(132) :: LineTmp,Line |

10 | |

11 |
INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na |

12 |
INTEGER(KINT) :: NbType, NbTypeUser |

13 |
INTEGER(KINT) :: I, J, Iat, NAtP, Idx |

14 |
INTEGER(KINT) :: IGeom |

15 |
REAL(KREAL) :: Xt,Yt,Zt |

16 |
LOGICAL :: Debug, TChk |

17 |
CHARACTER(4), ALLOCATABLE :: FFFt(:,:) ! 3,Na |

18 |
LOGICAL, ALLOCATABLE :: FFFl(:,:) ! 3,Na |

19 | |

20 |
!!! for VASP inputs |

21 |
! latx_loc used to check that the lattice parameters are the same for all geometries |

22 |
REAL(KREAL) :: lata_loc(3), latb_loc(3), latc_loc(3) |

23 |
! VaspPar_loc is the local value of the vasp parametre |

24 |
REAL(KREAL) :: VaspPar_loc |

25 |
! Vdir_loc is used to convert the read geometry to Cartesian. |

26 |
! V_direct is set by the first geometry. |

27 |
CHARACTER(LCHARS) :: Vdir_loc |

28 | |

29 | |

30 |
INTERFACE |

31 |
function valid(string) result (isValid) |

32 |
CHARACTER(*), intent(in) :: string |

33 |
logical :: isValid |

34 |
END function VALID |

35 |
END INTERFACE |

36 | |

37 |
debug=valid('Read_geom').or.valid('ReadGeom_vasp') |

38 | |

39 |
if (debug) Call Header("Entering ReadGeom_vasp") |

40 | |

41 |
ALLOCATE(FFF(3,nat),FFFl(3,Nat),FFFt(3,Nat)) |

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

43 |
! many things |

44 |
IGeom=1 |

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

46 |
READ(IOIN,'(A)') Vasp_Title |

47 |
READ(IOIN,*) Vasp_param |

48 | |

49 |
READ(IOIN,*) Lat_a |

50 |
READ(IOIN,*) Lat_b |

51 |
READ(IOIN,*) Lat_c |

52 | |

53 |
Lat_a=Lat_a*Vasp_param |

54 |
Lat_b=Lat_b*Vasp_param |

55 |
Lat_c=Lat_c*Vasp_param |

56 | |

57 |
ALLOCATE(NbAtType(nat)) |

58 |
READ(IOIN,'(A)') Vasp_types |

59 |
! First, how many different types ? |

60 |
NbAtType=0 |

61 |
READ(Vasp_types,*,END=999) NbAtType |

62 |
999 CONTINUE |

63 |
NbType=0 |

64 |
NatP=0 |

65 |
DO WHILE (NbAtType(NbType+1).NE.0) |

66 |
NbType=NbType+1 |

67 |
Natp=Natp+NbAtType(NbType) |

68 |
END DO |

69 | |

70 | |

71 |
if (NAtp.NE.Nat) THEN |

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

73 |
WRITE(IOOUT,*) "Using:",Natp |

74 |
DEALLOCATE(XyzGeomI, AtName) |

75 |
Nat=Natp |

76 |
ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt)) |

77 |
END if |

78 | |

79 |
! Do we know the atom types ? |

80 |
IF (AtTypes(1).EQ.' ') THEN |

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

82 |
! by looking into the POTCAR file... |

83 |
INQUIRE(File="POTCAR", EXIST=Tchk) |

84 |
IF (.NOT.Tchk) THEN |

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

86 |
STOP |

87 |
END IF |

88 |
OPEN(IOTMP,File="POTCAR") |

89 |
DO I=1,NbType |

90 |
Line='Empty' |

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

92 |
READ(IOTMP,'(A)') Line |

93 |
Line=AdjustL(Line) |

94 |
END DO |

95 |
Line=adjustl(Line(9:)) |

96 |
Idx=index(Line," ") |

97 |
Line=adjustl(Line(Idx+1:)) |

98 |
AtTypes(I)=Line(1:2) |

99 |
END DO |

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

101 |
CLOSE(IOTMP) |

102 | |

103 |
ELSE !AtTypes(1).EQ.' ' |

104 |
! user has provided atom types |

105 |
NbTypeUser=0 |

106 |
DO WHILE (AtTypes(NbTypeUser+1).NE.' ') |

107 |
NbTypeUser=NbTypeUser+1 |

108 |
END DO |

109 |
IF (NbType.NE.NbTypeUser) THEN |

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

111 |
STOP |

112 |
END IF |

113 |
END IF |

114 | |

115 |
IAt=1 |

116 |
DO I=1,NbType |

117 |
DO J=1,NbAtType(I) |

118 |
AtName(Iat)=AtTypes(I) |

119 |
Iat=Iat+1 |

120 |
END DO |

121 |
END DO |

122 |
DEALLOCATE(NbAtType) |

123 | |

124 |
NbTypes=NbType |

125 | |

126 |
READ(IOIN,'(A)') Vasp_comment |

127 |
READ(IOIN,'(A)') Vasp_direct |

128 | |

129 |
V_direct=adjustl(Vasp_direct) |

130 |
Call upcase(V_direct) |

131 | |

132 |
if (debug) WRITE(*,*) "V_direct=",Trim(V_Direct) |

133 | |

134 |
DO I=1,Nat |

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

136 |
DO J=1,3 |

137 |
FFF(J,I)=AdjustL(FFF(J,I)) |

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

139 |
END DO |

140 | |

141 |
END DO |

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

143 |
DO I=1,Nat |

144 |
Xt=XyzGeomI(IGeom,1,I) |

145 |
Yt=XyzGeomI(IGeom,2,I) |

146 |
Zt=XyzGeomI(IGeom,3,I) |

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

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

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

150 |
END DO |

151 | |

152 |
END IF |

153 | |

154 |
ALLOCATE(X0_vasp(Nat),Y0_vasp(Nat),Z0_vasp(Nat)) |

155 |
X0_vasp=XyzGeomI(IGeom,1,:) |

156 |
Y0_vasp=XyzGeomI(IGeom,2,:) |

157 |
Z0_vasp=XyzGeomI(IGeom,3,:) |

158 | |

159 |
if (debug) THEN |

160 |
WRITE(*,*) Nat |

161 |
WRITE(*,*) "Geom ",IGeom,"/",NGeomI |

162 |
DO I=1,Nat |

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

164 |
END DO |

165 |
END IF |

166 | |

167 | |

168 |
DO IGeom=2,NGeomI |

169 |
!!! PFL 16th March 2008 -> |

170 |
! Users might want to use different conditions for different geometries |

171 |
! i.e. to mix Cartesian and Direct :-p |

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

173 |
! |

174 | |

175 |
! Title |

176 |
READ(IOIN,*) LineTmp |

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

178 |
! Vasp_Param |

179 |
READ(IOIN,*) VaspPar_loc |

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

181 |
! Lattice parameters |

182 |
READ(IOIN,*) Lata_loc |

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

184 |
READ(IOIN,*) latb_loc |

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

186 |
READ(IOIN,*) Latc_loc |

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

188 |
! Attypes |

189 |
READ(IOIN,*) LineTmp |

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

191 |
! Comment |

192 |
READ(IOIN,*) LineTmp,line |

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

194 |
! V_direct |

195 |
READ(IOIN,*) Vdir_loc |

196 |
Vdir_loc=adjustl(Vdir_loc) |

197 |
Call upcase(Vdir_loc) |

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

199 | |

200 |
Lata_loc=Lata_loc*VaspPar_loc-Lat_a |

201 |
if (dot_product(lata_loc,Lata_loc).GE.1e-6) THEN |

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

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

204 |
STOP |

205 |
END IF |

206 | |

207 |
Latb_loc=Latb_loc*VaspPar_loc-Lat_b |

208 |
if (dot_product(latb_loc,Latb_loc).GE.1e-6) THEN |

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

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

211 |
STOP |

212 |
END IF |

213 | |

214 |
Latc_loc=Latc_loc*VaspPar_loc-Lat_c |

215 |
if (dot_product(latc_loc,Latc_loc).GE.1e-6) THEN |

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

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

218 |
STOP |

219 |
END IF |

220 | |

221 |
DO I=1,Nat |

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

223 |
DO J=1,3 |

224 |
! A frozen atom in any geom is considered frozen in all geom |

225 |
FFFt(j,I)=AdjustL(FFFt(j,i)) |

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

227 |
END DO |

228 |
END DO |

229 | |

230 |
IF (Vdir_loc(1:6).EQ.'DIRECT') THEN |

231 |
DO I=1,Nat |

232 |
Xt=XyzGeomI(IGeom,1,I) |

233 |
Yt=XyzGeomI(IGeom,2,I) |

234 |
Zt=XyzGeomI(IGeom,3,I) |

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

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

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

238 |
END DO |

239 |
END IF |

240 | |

241 |
if (debug) THEN |

242 |
WRITE(*,*) Nat |

243 |
WRITE(*,*) "Geom ",IGeom,"/",NGeomI |

244 |
DO I=1,Nat |

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

246 |
END DO |

247 |
END IF |

248 | |

249 |
END DO |

250 | |

251 |
DO I=1,Nat |

252 |
DO J=1,3 |

253 |
if (FFFl(j,i)) FFF(j,i)="T" |

254 |
END DO |

255 |
END DO |

256 | |

257 |
Deallocate(FFFl,FFFt) |

258 | |

259 | |

260 |
if (debug) Call Header("Exiting ReadGeom_vasp") |

261 | |

262 |
END SUBROUTINE ReadGeom_vasp |