Statistics
| Revision:

root / src / ReadInput_siesta.f90 @ 10

History | View | Annotate | Download (16.7 kB)

1
 SUBROUTINE ReadInput_siesta
2

    
3
! This routine reads an input template for Siesta
4

    
5
  use VarTypes
6
  use Path_module
7
  use Io_module
8

    
9
  IMPLICIT NONE
10

    
11
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
12
!
13
  INTERFACE
14
!
15
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16

    
17
     function valid(string) result (isValid)
18
       CHARACTER(*), intent(in) :: string
19
       logical                  :: isValid
20
     END function VALID
21

    
22
    FUNCTION SearchInput(Input,String,Line,Clean)  Result (Found)
23

    
24
      Use Vartypes
25
      use io_module
26
! Input
27
      TYPE (Input_line), POINTER, INTENT(IN) :: Input
28
      CHARACTER(*), INTENT(IN) :: String
29
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean
30

    
31
! Output
32
      TYPE (Input_line), POINTER, INTENT(OUT) :: Line
33

    
34
      LOGICAL :: Found
35

    
36
    END FUNCTION SearchInput
37

    
38
    FUNCTION InString(Line,String,Case,Clean,Back)  Result(Pos)
39

    
40
      Use VarTypes
41

    
42
      implicit none
43
! Input
44
      CHARACTER(*), INTENT(IN) :: Line
45
      CHARACTER(*), INTENT(IN) :: String
46
      LOGICAL, OPTIONAL, INTENT(IN) :: Case
47
      LOGICAL, OPTIONAL, INTENT(IN) :: Back
48
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean
49

    
50
! Output
51
! the position of String in Line (the first one) unless Back is present
52
      INTEGER(KINT) :: Pos
53
    END FUNCTION InString
54

    
55
    SUBROUTINE die(routine, msg, file, line, unit)
56

    
57
      Use VarTypes
58
      Use io_module
59

    
60
      implicit none
61
!--------------------------------------------------------------- Input Variables
62
      character(len=*), intent(in)           :: routine, msg
63
      character(len=*), intent(in), optional :: file
64
      integer(KINT), intent(in), optional      :: line, unit
65

    
66
    END SUBROUTINE die
67

    
68
    SUBROUTINE Warning(routine, msg, file, line, unit)
69

    
70
      Use VarTypes
71
      Use io_module
72

    
73
      implicit none
74

    
75
      character(len=*), intent(in)           :: routine, msg
76
      character(len=*), intent(in), optional :: file
77
      integer(KINT), intent(in), optional      :: line, unit
78

    
79
    END SUBROUTINE Warning
80

    
81

    
82
 SUBROUTINE WriteList(Input,Unit)
83

    
84
! This routine reads an input template for Siesta
85

    
86
  use VarTypes
87
  use Io_module
88

    
89
  IMPLICIT NONE
90

    
91
! Input
92
      TYPE (Input_line), POINTER, INTENT(IN) :: Input
93
      INTEGER(KINT), OPTIONAL, INTENT(IN) :: Unit
94

    
95
    END SUBROUTINE WriteList
96

    
97
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
98
!
99
 END  INTERFACE
100
!
101
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
102

    
103

    
104
  CHARACTER(132) ::  Line,Line2
105
  INTEGER(KINT) ::  Idx,LineL
106
  INTEGER(KINT) ::  IAt,I
107
  INTEGER(KINT) :: IoRead, ITmp
108
  REAL(KREAL) :: B(3),Xtmp, Ytmp, Ztmp
109
  REAL(KREAL) :: Lata,Latb,Latc,Alpha,Beta,Gamma
110
  REAL(KREAL), ALLOCATABLE :: GeomCartTmp(:,:) ! (Nat,3)
111
  REAL(KREAL) :: LatLoc(3,3)
112
  
113
  LOGICAL :: Debug
114
  LOGICAL :: FSpecies, FCoord
115

    
116
  TYPE(Input_Line), POINTER :: Search,Bla
117

    
118
  Debug=Valid("readinput").OR.Valid("readinput_siesta")
119

    
120
 if (debug) Call Header("Entering ReadInput_Siesta")
121

    
122
  ! We read the Siesta input file
123

    
124
     ALLOCATE(Siesta_Input)
125
     NULLIFY(Siesta_Input%next)
126
     NULLIFY(Siesta_Input%prev)
127
     Current => Siesta_Input
128

    
129
     FSpecies=.FALSE.
130
     FCoord=.FALSE.
131

    
132
     READ(IOIN,'(A)',iostat=IoRead) Line
133

    
134
     DO WHILE (IoRead==0)
135
        Line=AdjustL(Line)
136
        if (debug) WRITE(*,*) 'Line:', Line
137
        current%Line=TRIM(Line)
138
        ALLOCATE(current%next)
139
        Current%next%prev => Current
140
        Current => Current%next
141
        current%line="toto"
142
        Nullify(Current%next)
143

    
144
        READ(IOIN,'(A)',iostat=IoRead) Line
145
     END DO
146
! With this procedure, the last current is not valid and should be deleted.
147
       Bla => Current%prev
148
       nullify(Bla%next)
149
       deallocate(Current)
150
       current => Bla
151

    
152
     if (debug) THEN
153
        WRITE(*,*) 'Input just read'
154
        Call WriteList(Siesta_input,Unit=IOOUT)
155
     END IF
156

    
157
! We analyse the input
158

    
159
! We look for the NumberofAtoms
160
     If (SearchInput(Siesta_Input,"NUMBEROFATOMS",Search,Clean=".-_")) THEN
161
        Line=AdjustL(Search%line)
162
        Idx=Index(Line," ")
163
        Line2=Trim(AdjustL(Line(Idx+1:)))
164
        Read(Line2,*) IAt
165
        if (Iat/=Nat) THEN
166
           Call Die('ReadInput_siesta','Nat in FDF sample different from  Nat Path input', Unit=IOOUT)
167
        END IF
168
     ELSE
169
! There is no atom number defined !!!
170
           Call Warning('ReadInput_siesta','No NumberOfAtoms in FDF sample', Unit=IOOUT)
171
           WRITE(current%Line,'(1X,"NumberOfAtoms ",I5)') Nat
172
           ALLOCATE(current%next)
173
           Current%next%prev => Current
174
           Current => Current%next
175
           Nullify(Current%next)
176
         
177
        END IF
178

    
179
! We look for the NumberofSpecies
180
        IF (SearchInput(Siesta_Input,"NUMBEROFSPECIES",Search,Clean=".-_")) THEN
181
           Line=AdjustL(Search%line)
182
           Idx=Index(Line," ")
183
           Line2=Trim(AdjustL(Line(Idx+1:)))
184
           Read(Line2,*) Siesta_NbSpecies
185
           ALLOCATE(Siesta_SpeciesMass(Siesta_NbSpecies))
186
           ALLOCATE(Siesta_SpeciesName(Siesta_NbSpecies))
187
        END IF
188

    
189
! We look for SystemLabel
190
        If (SearchInput(Siesta_Input,"SYSTEMLABEL",Search,Clean=".-_")) THEN
191
           Line=AdjustL(Search%line)
192
           Idx=Index(Line," ")
193
           Siesta_Label=Trim(adjustl(Line(Idx+1:)))
194
        ELSE
195
           Siesta_label='siesta'
196
        END IF
197

    
198
! We look for the ChemicalSpeciesLabel block
199
        IF (SearchInput(Siesta_Input,"CHEMICALSPECIESLABEL",Search,Clean=".-_")) THEN
200
           if (.NOT.ALLOCATED(IdxSpecies)) ALLOCATE(IdxSpecies(NAt))
201
           I=0
202
           Search => Search%next
203
           DO WHILE (InString(Search%line,'ENDBLOCK',Case=.FALSE.,Clean=".-_")==0)
204
              I=I+1
205
              if (I>Siesta_NbSpecies) THEN
206
                 Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel','Found more line in this block than NbSpecies !')
207
              END IF
208
              Line=AdjustL(Search%Line)
209
              Read(Line,*) ITmp
210
              Idx=Index(Line,' ')
211
              Line=AdjustL(Line(Idx+1:))
212
              Read(Line,*) Idx
213
              Siesta_SpeciesMass(ITmp)=Idx
214
              Idx=Index(Line,' ')
215
              Line=AdjustL(Line(Idx+1:))
216
              Idx=Index(Line,' ')
217
              Siesta_SpeciesName(ITmp)=Line(1:Idx-1)
218
              Search => Search%next
219
           END DO
220
           IF (I/=Siesta_NbSpecies) Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel', &
221
                'Number of lines in this block different from NbSpecies')
222
        ELSE
223
           Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel', &
224
                'Block ChemicalSpeciesLabel is mandatory !')
225
        END IF
226

    
227

    
228
! We look for the AtomicCoordinatesAndAtomicSpecies  block
229
        IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESANDATOMICSPECIES",Search,Clean=".-_")) THEN
230
           ALLOCATE(Siesta_Paste(Nat))
231
           Current=>Search
232
           DO I=1,NAt
233
              Search => Search%next
234
              Read(Search%line,*) XTmp,YTmp,ZTmp,IdxSpecies(I)
235
! We give a name to this atom
236
              IF (AtName(I)/="") THEN
237
!                    Write(Line,'(A,I5,1X,A2,1X,A32)') 'AtName & SpeciesName for atom ',I,AtName(I), &
238
!                         TRIM(Siesta_SpeciesName(IdxSpecies(I)))
239
!                    WRITE(*,'(1X,A)') Line
240
                 If (InString(Atname(I),TRIM(Siesta_SpeciesName(IdxSpecies(I))))==0) THEN
241
                    Write(Line,'(A,I5,1X,A2,1X,A32)') 'AtName /= SpeciesName for atom ',I,AtName(I), &
242
                         TRIM(Siesta_SpeciesName(IdxSpecies(I)))
243
                    Call Die('Readinput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)
244
                 END IF
245
! We look for the Mass number also
246
                 IF (Atome(I)/=Siesta_SpeciesMass(IdxSpecies(I))) THEN
247
                    Write(Line,'(A,I5,1X,I5,1X,I5)') 'AtMass /= SpeciesMass for atom ',I,Atome(I), &
248
                         Siesta_SpeciesMass(IdxSpecies(I))
249
                    Call Die('ReadInput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)
250
                 END IF
251
              ELSE
252
                 AtName(I)=AdjustL(Trim(Siesta_SpeciesName(IdxSpecies(I))))
253
                 Atome(I)=Siesta_SpeciesMass(IdxSpecies(I))
254
              END IF
255
! we look for something else at the end of the line
256
! We save everything but the x,y,z, and species description              
257
              Line=AdjustL(search%line)
258
! We skip x
259
              Idx=Index(Line,' ')
260
              Line=AdjustL(Line(Idx+1:))
261
! we skip y
262
              Idx=Index(Line,' ')
263
              Line=AdjustL(Line(Idx+1:))
264
! we skip z
265
              Idx=Index(Line,' ')
266
              Line=AdjustL(Line(Idx+1:))
267
! we skip species
268
              Idx=Index(Line,' ')
269
              Line=AdjustL(Line(Idx+1:))
270
              Siesta_Paste(I)=TRIM(Line)
271
           END DO
272
           if (debug) THEN
273
              WRITE(*,*) 'Input before deletion of coord block'
274
              Call WriteList(Siesta_input)
275
           END IF
276
! We will now delete this block from our input sample as it will then be
277
! written directly by Opt'n Path
278
! Search%next point on %endblock
279
           search => search%next
280
           IF (ASSOCIATED(Search%next)) THEN
281
! we are not at the end of the input file
282
              if (ASSOCIATED(Current%Prev)) THEN
283
                 Search%next%prev => current%prev
284
                 current%prev%next => search%next
285
              ELSE
286
! the coordinate block is at the begining of the input
287
                 siesta_input => Search%next
288
                 Nullify(Siesta_Input%Prev)
289
              END IF
290
           ELSE
291
! the coordinate block is the last part of the input
292
              nullify(current%prev%next)
293
           END IF
294
           if (debug) THEN
295
              WRITE(*,*) 'Input after deletion of coord block 1'
296
              Call WriteList(Siesta_input)
297
           END IF
298

    
299
           DO I=1,Nat+2
300
              bla=>current
301
              current=> current%next
302
              deallocate(bla)
303
           END DO
304
           if (debug) THEN
305
              WRITE(*,*) 'Input after deletion of coord block 2'
306
              Call WriteList(Siesta_input)
307
           END IF
308
        ELSE
309
           IF (SearchInput(Siesta_Input,"ZMATRIX",Search,Clean=".-_")) THEN
310
              call Die('ReadInput_Siesta','For now, I need the full block'//&
311
                   'AtomicCoordinatesAndAtomicSpecies, ZMatrix not yet handled.')
312
           ELSE
313
              call Die('ReadInput_Siesta','For now, I need the full block' //&
314
              'AtomicCoordinatesAndAtomicSpecies.')
315
           END IF
316
        END IF
317

    
318
        IF (SearchInput(Siesta_Input,"LATTICE",Search,Clean=".-_")) THEN
319
! This is a pbc calculation
320
           FPBC=.TRUE.
321
           IPer=3
322
           Kabeg=-1
323
           kaEnd=1
324
           kbBeg=-1
325
           kbEnd=1
326
           kcBeg=-1
327
           kcEnd=1
328
         
329

    
330
         IF (SearchInput(Siesta_Input,"LATTICECONSTANT",Search,Clean=".-_")) THEN
331
           Line=Adjustl(Search%Line)
332
! We discar the label
333
           Idx=Index(Line,' ')
334
           Line=AdjustL(Line(Idx+1:))
335
! we read the value
336
           Read(Line,*) Siesta_LatticeConstant
337
! We discard the value
338
           Idx=Index(Line,' ')
339
           Line=AdjustL(Line(Idx+1:))
340
           LineL=Len_Trim(Line)
341
! We read the unit
342
           If (LineL>0) THEN
343
              Call UpCase(Line)
344
              SELECT CASE (Line) 
345
                CASE ('ANG')
346
                   Siesta_Lat_Unit=1.d0
347
                CASE ('BOHR')
348
                   Siesta_Lat_Unit=a0
349
              END SELECT
350
           ELSE
351
              Siesta_Lat_Unit=1.d0
352
           END IF
353
         ELSE
354
! We have no LatticeConstant but we need one for our calculation
355
            Siesta_LatticeConstant=1.
356
            Siesta_Lat_Unit=1.d0
357
         END IF
358

    
359
         IF (SearchInput(Siesta_Input,"LATTICEVECTORS",Search,Clean=".-_")) THEN
360
           search => search%next
361
           Line=Adjustl(Search%Line)
362
           Read(Line,*) Lat_a(1:3)
363
           search => search%next
364
           Line=Adjustl(Search%Line)
365
           Read(Line,*) Lat_b(1:3)
366
           search => search%next
367
           Line=Adjustl(Search%Line)
368
           Read(Line,*) Lat_c(1:3)
369

    
370
         ELSEIF (SearchInput(Siesta_Input,"LATTICEPARAMETERS",Search,Clean=".-_")) THEN
371
           search => search%next
372
           Line=Adjustl(Search%Line)
373
           Read(Line,*) Lata,Latb,Latc,Alpha,Beta,Gamma
374

    
375
            Lat_a=0.
376
            Lat_b=0.
377
            Lat_c=0.
378
            Lat_a(1)=Lata
379
            Lat_b(1)=LatB*cos(Gamma*Pi/180d0)
380
            Lat_b(2)=LatB*sin(Gamma*Pi/180d0)
381
            Lat_c(1)=Latc*cos(Beta*Pi/180.)
382
            Lat_c(2)=(Latb*Latc*cos(Alpha*Pi/180.)-Lat_b(1)*Lat_c(1))/Lat_b(2)
383
            Lat_c(3)=sqrt(Latc*Latc-Lat_c(1)**2-Lat_c(2)**2)
384
         ELSE
385
! We have a lattice constant but nothing else. We issue a warning just in case,
386
! and we take the defaut values 1. 1. 1. 90. 90. 90.
387
            Line="LatticeConstant found, but not LatticeVectors nor LatticeParameters. Taking 1. 1. 1. 90. 90. 90."
388
            Call Warning("ReadInput_siesta",Line)
389
            Lat_a=(/1.,0.,0./)
390
            Lat_b=(/0.,1.,0./)
391
            Lat_c=(/0.,0.,1./)
392
         END IF
393

    
394
         IF (debug) THEN 
395
            WRITE(*,*) "Lattice vectors before scaling with LatticeConstant"
396
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_a",Lat_a
397
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_b",Lat_b
398
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_c",Lat_c
399
         END IF
400
           Lat_a=Lat_a*Siesta_LatticeConstant*Siesta_Lat_Unit
401
           Lat_b=Lat_b*Siesta_LatticeConstant*Siesta_Lat_Unit
402
           Lat_c=Lat_c*Siesta_LatticeConstant*Siesta_Lat_Unit
403

    
404
         IF (debug) THEN 
405
            WRITE(*,*) "Lattice vectors After scaling with LatticeConstant"
406
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_a",Lat_a
407
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_b",Lat_b
408
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_c",Lat_c
409
         END IF
410

    
411
      END IF
412

    
413
         IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESFORMAT",Search,Clean=".-_")) THEN
414
           Line=Adjustl(Search%Line)
415
           Idx=Index(Line,' ')
416
           Line=AdjustL(Line(Idx+1:))
417
           Call UpCase(Line)
418
           SELECT CASE (Line) 
419
              CASE ('ANG','NOTSCALEDCARTESIANANG')
420
                 Siesta_Unit_Read=1.d0
421
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
422
                 IF (FPBC) THEN
423
                    Siesta_Unit_Read=1.d0
424
                    V_DIRECT="DIRECT"
425
                    Latr(1:3,1)=Lat_a
426
                    Latr(1:3,2)=Lat_b
427
                    Latr(1:3,3)=Lat_c
428
                    B=1.
429
! TO DO: replace by LAPACK
430
                    CALL Gaussj(Latr,3,3,B,1,1)
431
                    LatLoc(1:3,1)=Lat_a
432
                    LatLoc(1:3,2)=Lat_b
433
                    LatLoc(1:3,3)=Lat_c
434
                    Allocate(GeomCartTmp(Nat,3))
435
                    Do I=1,NGeomI
436
                       GeomCartTmp=Reshape(XyzGeomI(I,1:3,1:Nat),(/Nat,3/),ORDER=(/2,1/))
437
! TO DO: shall we replace MatMull by DGEMM ?
438
                       XyzGeomI(I,1:3,1:Nat)=Reshape(MatMul(GeomCartTmp,Latloc),(/3,Nat/),ORDER=(/2,1/))
439
                   END DO
440
                 ELSE
441
                    WRITE(Line2,*) "AtomicCoordinatesFormat=" // Trim(Line) &
442
               //", but there is no information about Lattice parameters"
443
                    Call Die("ReadInput_siesta",Line2)
444
                 END IF
445
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
446
                 Siesta_Unit_Read=a0
447
                 IF (INPUT=='SIESTA') THEN
448
! We have read the coordinates, but not the unit. This is corrected here
449
                    XyZGeomI=XyzGeomI*a0
450
                 END IF
451
              CASE ('SCALEDCARTESIAN')
452
                 Siesta_Unit_Read=Siesta_LatticeConstant
453
                 If (INPUT=='SIESTA') THEN
454
                    XyzGeomI=XyzGeomI*Siesta_Unit_Read
455
                 END IF
456
            END SELECT
457
         ELSE
458
            Siesta_Unit_Read=1.d0
459
         END IF
460

    
461
        IF (SearchInput(Siesta_Input,"ATOMICCOORFORMATOUT",Search,Clean=".-_")) THEN
462
           Line=Adjustl(Search%Line)
463
           Idx=Index(Line,' ')
464
           Line=AdjustL(Line(Idx+1:))
465
           Call UpCase(Line)
466
           SELECT CASE (Line) 
467
              CASE ('ANG','NOTSCALEDCARTESIANANG')
468
                 Siesta_Unit_Write=1.d0
469
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
470
                 Siesta_Unit_Write=1.d0
471
                 V_DIRECT_Write='DIRECT'
472
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
473
                 Siesta_Unit_Write=Unit ! 1/a0
474
            END SELECT
475
         END IF
476
         IF (SearchInput(Siesta_Input,"SUPERCELL",Search,Clean=".-_")) THEN
477
            Line=Adjustl(search%line)
478
            Call Upcase(Line)
479
            Call CleanString(Line,".-_")
480
            If (Index(Line,"%BLOCK")/=0) THEN
481
               Call Die('ReadInput_siesta:SuperCell',"SuperCell  NOT possible in Opt'n Path")
482
            END IF
483
         END IF
484

    
485
!         Call Die('Readinput_siesta:fin',' Lecture finie')
486

    
487
     if (debug) Call Header("Exiting ReadInput_Siesta")
488

    
489
   END SUBROUTINE READINPUT_SIESTA