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