Statistiques
| Révision :

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