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 $user-defined, 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_loc-Lat_a
255
        if (dot_product(lata_loc,Lata_loc).GE.1e-6) 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_loc-Lat_b
262
        if (dot_product(latb_loc,Latb_loc).GE.1e-6) 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_loc-Lat_c
269
        if (dot_product(latc_loc,Latc_loc).GE.1e-6) 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