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*LatcLat_c(1)**2Lat_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