Statistics
| Revision:

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