Revision 10 src/ReadInput_siesta.f90

ReadInput_siesta.f90 (revision 10)
102 102

  
103 103

  
104 104
  CHARACTER(132) ::  Line,Line2
105
  INTEGER(KINT) ::  Idx
105
  INTEGER(KINT) ::  Idx,LineL
106 106
  INTEGER(KINT) ::  IAt,I
107 107
  INTEGER(KINT) :: IoRead, ITmp
108
  REAL(KREAL) :: Xtmp, Ytmp, Ztmp
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)
109 112
  
110 113
  LOGICAL :: Debug
111 114
  LOGICAL :: FSpecies, FCoord
......
312 315
           END IF
313 316
        END IF
314 317

  
315
        IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESFORMAT",Search,Clean=".-_")) THEN
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
316 331
           Line=Adjustl(Search%Line)
332
! We discar the label
317 333
           Idx=Index(Line,' ')
318 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:))
319 417
           Call UpCase(Line)
320 418
           SELECT CASE (Line) 
321 419
              CASE ('ANG','NOTSCALEDCARTESIANANG')
322 420
                 Siesta_Unit_Read=1.d0
323 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
324 445
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
325 446
                 Siesta_Unit_Read=a0
326 447
                 IF (INPUT=='SIESTA') THEN
327 448
! We have read the coordinates, but not the unit. This is corrected here
328 449
                    XyZGeomI=XyzGeomI*a0
329 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
330 456
            END SELECT
457
         ELSE
458
            Siesta_Unit_Read=1.d0
331 459
         END IF
332 460

  
333 461
        IF (SearchInput(Siesta_Input,"ATOMICCOORFORMATOUT",Search,Clean=".-_")) THEN
......
339 467
              CASE ('ANG','NOTSCALEDCARTESIANANG')
340 468
                 Siesta_Unit_Write=1.d0
341 469
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
470
                 Siesta_Unit_Write=1.d0
471
                 V_DIRECT_Write='DIRECT'
342 472
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
343
                 Siesta_Unit_Write=Unit
473
                 Siesta_Unit_Write=Unit ! 1/a0
344 474
            END SELECT
345 475
         END IF
346

  
347
         IF (SearchInput(Siesta_Input,"LATTICECONSTANT",Search,Clean=".-_")) THEN
348
           Line=Adjustl(Search%Line)
349
! We discar the label
350
           Idx=Index(Line,' ')
351
           Line=AdjustL(Line(Idx+1:))
352
! we read the value
353
           Read(Line,*) Siesta_LatticeConstant
354
! We discard the value
355
           Idx=Index(Line,' ')
356
           Line=AdjustL(Line(Idx+1:))
357
! We read the unit
358
           Call UpCase(Line)
359
           SELECT CASE (Line) 
360
              CASE ('ANG')
361
                 Siesta_Lat_Unit=1.d0
362
              CASE ('BOHR')
363
                 Siesta_Lat_Unit=Unit
364
            END SELECT
365
! for now!
366
            Call Die('ReadInput_siesta:LatticeConstant',"For now, periodic calculations are NOT possible in Opt'n Path")
367

  
368
         END IF
369

  
370
         IF (SearchInput(Siesta_Input,"LATTICEVECTORS",Search,Clean=".-_")) THEN
371

  
372
            Call Die('ReadInput_siesta:LatticeVectors',"For now, periodic calculations are NOT possible in Opt'n Path")
373

  
374
         END IF
375

  
376
         IF (SearchInput(Siesta_Input,"LATTICEPARAMETERS",Search,Clean=".-_")) THEN
377

  
378
            Call Die('ReadInput_siesta:LatticeParameters',"For now, periodic calculations are NOT possible in Opt'n Path")
379

  
380
         END IF
381

  
382 476
         IF (SearchInput(Siesta_Input,"SUPERCELL",Search,Clean=".-_")) THEN
383

  
384
            Call Die('ReadInput_siesta:SuperCell',"SuperCell  NOT possible in Opt'n Path")
385

  
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
386 483
         END IF
387 484

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

Also available in: Unified diff