Révision 10
| src/egrad_gaussian.f90 (revision 10) | ||
|---|---|---|
| 93 | 93 |
DO I=1,Nat |
| 94 | 94 |
If (renum) THEN |
| 95 | 95 |
Iat=Order(I) |
| 96 |
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),GeomCart(Iat,:), TRIM(Gauss_Paste(I)) |
|
| 96 |
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),1X,A)') Trim(AtName(I)),GeomCart(Iat,:), TRIM(Gauss_Paste(I))
|
|
| 97 | 97 |
ELSE |
| 98 | 98 |
Iat=OrderInv(I) |
| 99 |
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomCart(I,:), TRIM(Gauss_Paste(Iat)) |
|
| 99 |
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),1X,A)') Trim(AtName(Iat)),GeomCart(I,:), TRIM(Gauss_Paste(Iat))
|
|
| 100 | 100 |
END IF |
| 101 | 101 |
END DO |
| 102 | 102 |
|
| src/egrad_mopac.f90 (revision 10) | ||
|---|---|---|
| 3 | 3 |
! This routines calculates the energy and the gradient of |
| 4 | 4 |
! a molecule, using mopac |
| 5 | 5 |
|
| 6 |
use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, ProgExe |
|
| 6 |
use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, ProgExe, & |
|
| 7 |
FPBC, Lat_a, Lat_b, Lat_c, IPer |
|
| 7 | 8 |
use Io_module |
| 8 | 9 |
|
| 9 | 10 |
IMPLICIT NONE |
| ... | ... | |
| 92 | 93 |
WRITE(IOTMP,'(1X,A10,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:) |
| 93 | 94 |
END IF |
| 94 | 95 |
END DO |
| 96 |
! If we have periodic boundaries, we have to print them |
|
| 97 |
IF (FPBC) THEN |
|
| 98 |
IF (IPER>=1) THEN |
|
| 99 |
WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_a(1:3) |
|
| 100 |
END IF |
|
| 101 |
IF (IPER>=2) THEN |
|
| 102 |
WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_b(1:3) |
|
| 103 |
END IF |
|
| 104 |
IF (IPER>=3) THEN |
|
| 105 |
WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_c(1:3) |
|
| 106 |
END IF |
|
| 107 |
END IF |
|
| 95 | 108 |
|
| 96 | 109 |
WRITE(IOTMP,'(1X,A)') Trim(Mopac_EndGeom) |
| 97 | 110 |
|
| src/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')
|
| src/ReadInput_mopac.f90 (revision 10) | ||
|---|---|---|
| 13 | 13 |
CHARACTER(*), intent(in) :: string |
| 14 | 14 |
logical :: isValid |
| 15 | 15 |
END function VALID |
| 16 |
|
|
| 17 |
SUBROUTINE die(routine, msg, file, line, unit) |
|
| 18 |
|
|
| 19 |
Use VarTypes |
|
| 20 |
Use io_module |
|
| 21 |
|
|
| 22 |
implicit none |
|
| 23 |
character(len=*), intent(in) :: routine, msg |
|
| 24 |
character(len=*), intent(in), optional :: file |
|
| 25 |
integer(KINT), intent(in), optional :: line, unit |
|
| 26 |
|
|
| 27 |
END SUBROUTINE die |
|
| 28 |
|
|
| 29 |
|
|
| 16 | 30 |
END INTERFACE |
| 17 | 31 |
|
| 18 | 32 |
|
| 19 |
CHARACTER(132) :: Line,Line2
|
|
| 33 |
CHARACTER(LCHARS) :: Line,LineUp
|
|
| 20 | 34 |
INTEGER(KINT) :: LineL, Idx, NTmp |
| 35 |
INTEGER(KINT) :: NatMopac |
|
| 36 |
REAL(KREAL) :: Lat(3,3) |
|
| 21 | 37 |
|
| 22 | 38 |
LOGICAL :: Debug |
| 23 | 39 |
|
| ... | ... | |
| 50 | 66 |
LineL=len(Trim(Line)) |
| 51 | 67 |
IF (Line(1:1)/="*") THEN |
| 52 | 68 |
IF (NTmp==0) THEN |
| 53 |
Line2=Line
|
|
| 54 |
Call UpCase(Line2)
|
|
| 55 |
Idx=Index(Line2,'GRADIENTS')
|
|
| 69 |
LineUp=Line
|
|
| 70 |
Call UpCase(LineUp)
|
|
| 71 |
Idx=Index(LineUp,'GRADIENTS')
|
|
| 56 | 72 |
If (Idx==0) Line=TRIM(Line) // " GRADIENTS" |
| 57 |
Idx=Index(Line2,'1SCF')
|
|
| 73 |
Idx=Index(LineUp,'1SCF')
|
|
| 58 | 74 |
If (Idx==0) Line=TRIM(Line) // " 1SCF" |
| 59 | 75 |
END IF |
| 60 | 76 |
current%Line=TRIM(Line) |
| ... | ... | |
| 77 | 93 |
! END DO |
| 78 | 94 |
|
| 79 | 95 |
! Now the geometry... that we just skip :) |
| 96 |
! PFL 2013 Apr |
|
| 97 |
! We take care that there is no Translation vectors... |
|
| 98 |
! We also check that the number of atoms is ok |
|
| 80 | 99 |
IF (DEBUG) WRITE(*,*) "Reading Mopac Geometry" |
| 81 | 100 |
Mopac_EndGeom="" |
| 82 | 101 |
LineL=1 |
| 102 |
NatMopac=0 |
|
| 103 |
Lat=0.d0 |
|
| 104 |
IPer=0 |
|
| 105 |
FPBC=.FALSE. |
|
| 83 | 106 |
DO WHILE (LineL.NE.0) |
| 84 | 107 |
READ(IOIN,'(A)',END=989) Line |
| 85 | 108 |
Line=AdjustL(Line) |
| 86 | 109 |
LineL=len(Trim(Line)) |
| 87 | 110 |
! The last line might be either blank or filled with 0 |
| 88 |
IF (Line(1:1)=="0") THEN |
|
| 89 |
LineL=0 |
|
| 90 |
Mopac_EndGeom=Trim(Line) |
|
| 111 |
If (LineL>0) THEN |
|
| 112 |
SELECT CASE (Line(1:1)) |
|
| 113 |
CASE ("0")
|
|
| 114 |
LineL=0 |
|
| 115 |
Mopac_EndGeom=Trim(Line) |
|
| 116 |
CASE("*")
|
|
| 117 |
CurCom%Line=TRIM(LINE) |
|
| 118 |
ALLOCATE(CurCom%Next) |
|
| 119 |
CurCom => CurCom%Next |
|
| 120 |
NULLIFY(CurCom%Next) |
|
| 121 |
CASE DEFAULT |
|
| 122 |
LineUp=Line |
|
| 123 |
Call UpCase(LineUp) |
|
| 124 |
If (LineUp(1:2)=="TV") THEN |
|
| 125 |
FPBC=.TRUE. |
|
| 126 |
IPer=IPer+1 |
|
| 127 |
If (Iper>3) THEN |
|
| 128 |
Call Die("ReadInput Mopac","Iper>3",Unit=IOOUT)
|
|
| 129 |
END IF |
|
| 130 |
NTmp=Index(LineUp," ") |
|
| 131 |
LineUp=LineUp(NTmp:) |
|
| 132 |
Read(LineUp,*) Lat(IPer,1:3) |
|
| 133 |
ELSE |
|
| 134 |
NatMopac=NatMopac+1 |
|
| 135 |
END IF |
|
| 136 |
END SELECT |
|
| 137 |
END IF |
|
| 138 |
|
|
| 139 |
END DO |
|
| 140 |
|
|
| 141 |
! WRITE(*,*) "NatMopac,Nat:",NAtMopac,Nat |
|
| 142 |
IF (NatMopac/=Nat) Call Die("ReadInput_mopac","Nat read does not mat nat",Unit=IOOUT)
|
|
| 143 |
IF (FPBC) THEN |
|
| 144 |
Lat_a(1:3)=Lat(1,1:3) |
|
| 145 |
Lat_b(1:3)=Lat(2,1:3) |
|
| 146 |
Lat_c(1:3)=Lat(3,1:3) |
|
| 147 |
If (IPer>=1) THEN |
|
| 148 |
kaBeg=-1 |
|
| 149 |
kaEnd=1 |
|
| 91 | 150 |
END IF |
| 92 |
IF (Line(1:1)=="*") THEN |
|
| 93 |
CurCom%Line=TRIM(LINE) |
|
| 94 |
ALLOCATE(CurCom%Next) |
|
| 95 |
CurCom => CurCom%Next |
|
| 96 |
NULLIFY(CurCom%Next) |
|
| 151 |
If (IPer>=2) THEN |
|
| 152 |
kbBeg=-1 |
|
| 153 |
kbEnd=1 |
|
| 97 | 154 |
END IF |
| 98 |
END DO |
|
| 155 |
If (IPer==3) THEN |
|
| 156 |
kcBeg=-1 |
|
| 157 |
kcEnd=1 |
|
| 158 |
END IF |
|
| 159 |
If (IPer>3) THEN |
|
| 160 |
Call Die("Readinput_mopac","Found too many Tv lines !",Unit=IOOUT)
|
|
| 161 |
END IF |
|
| 162 |
END IF |
|
| 99 | 163 |
|
| 100 | 164 |
! If we are here, there might be something else to read: Mopac_end |
| 101 | 165 |
|
| src/Io_module.f90 (revision 10) | ||
|---|---|---|
| 10 | 10 |
INTEGER(KINT) :: IOIN=11, IOOUT=12, IOCART=14 |
| 11 | 11 |
INTEGER(KINT) :: IOGEOM=15, IODAT=16,IoGplot=17 |
| 12 | 12 |
INTEGER(KINT), PARAMETER :: IOTMP=21,IOTMP2=22, IOTMP3=23 |
| 13 |
INTEGER(KINT), PARAMETER :: IOERR=19
|
|
| 13 |
INTEGER(KINT), PARAMETER :: IOERR=0
|
|
| 14 | 14 |
CHARACTER(SCHARS) :: RunMode |
| 15 | 15 |
|
| 16 | 16 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| src/ReadInput_vasp.f90 (revision 10) | ||
|---|---|---|
| 33 | 33 |
|
| 34 | 34 |
if (debug) Call Header("Entering ReadInput_Vasp")
|
| 35 | 35 |
|
| 36 |
FPBC=.TRUE. |
|
| 37 |
|
|
| 36 | 38 |
if (Input/="VASP") THEN |
| 37 | 39 |
|
| 38 | 40 |
! Input was not Vasp, so many parameters are missing like lattice |
| ... | ... | |
| 138 | 140 |
Latr(1:3,2)=Lat_b |
| 139 | 141 |
Latr(1:3,3)=Lat_c |
| 140 | 142 |
B=1. |
| 143 |
! TO DO: replace by LAPACK |
|
| 141 | 144 |
CALL Gaussj(Latr,3,3,B,1,1) |
| 142 | 145 |
ELSE |
| 143 | 146 |
Latr=0. |
| ... | ... | |
| 180 | 183 |
Vasp_Title=Trim(Line) // " " // Trim(adjustL(Vasp_Title)) |
| 181 | 184 |
if (debug) WRITE(*,*) "DBG MAIN, L809, VMD=T, Vasp_Title=",TRIM(Vasp_Title) |
| 182 | 185 |
END IF |
| 183 |
Call CheckPeriodicBound |
|
| 184 | 186 |
|
| 185 | 187 |
if (debug) Call Header("Exiting ReadInput_Vasp")
|
| 186 | 188 |
|
| src/Path.f90 (revision 10) | ||
|---|---|---|
| 280 | 280 |
|
| 281 | 281 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| 282 | 282 |
|
| 283 |
LOGICAL :: Debug, Fini,Flag_Opt_Geom |
|
| 283 |
LOGICAL :: Debug, Fini,Flag_Opt_Geom,TChk
|
|
| 284 | 284 |
|
| 285 | 285 |
! INTEGER(KINT), EXTERNAL :: Iargc |
| 286 | 286 |
INTEGER(KINT), EXTERNAL :: ConvertNumAt |
| ... | ... | |
| 831 | 831 |
If (ProgExe=="EMPTY") ProgExe="VASP" |
| 832 | 832 |
UnitProg="eV" |
| 833 | 833 |
ConvE=eV2kcal |
| 834 |
FPBC=.TRUE. |
|
| 835 |
IPer=3 |
|
| 836 |
Kabeg=-1 |
|
| 837 |
kaEnd=1 |
|
| 838 |
kbBeg=-1 |
|
| 839 |
kbEnd=1 |
|
| 840 |
kcBeg=-1 |
|
| 841 |
kcEnd=1 |
|
| 834 | 842 |
CASE ("TM","TURBOMOLE")
|
| 835 | 843 |
Prog="TURBOMOLE" |
| 836 | 844 |
If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep" |
| 837 | 845 |
UnitProg="au" |
| 838 | 846 |
ConvE=au2kcal |
| 847 |
FPBC=.FALSE. |
|
| 839 | 848 |
CASE ("TEST","CHAMFRE","TEST2D","LEPS")
|
| 840 | 849 |
ProgExe="" |
| 841 | 850 |
UnitProg="au" |
| ... | ... | |
| 1192 | 1201 |
|
| 1193 | 1202 |
Call ReadInput |
| 1194 | 1203 |
|
| 1204 |
if (FPBC) THEN |
|
| 1205 |
Allocate(XGeomRefPBC(Nat),YGeomRefPBC(Nat),ZGeomRefPBC(Nat)) |
|
| 1206 |
XGeomRefPBC(1:Nat)=XyzGeomI(1,1,1:Nat) |
|
| 1207 |
YGeomRefPBC(1:Nat)=XyzGeomI(1,2,1:Nat) |
|
| 1208 |
ZGeomRefPBC(1:Nat)=XyzGeomI(1,3,1:Nat) |
|
| 1209 |
|
|
| 1210 |
Call CheckPeriodicBound |
|
| 1211 |
END IF |
|
| 1212 |
|
|
| 1195 | 1213 |
IF (COORD=="MIXED") Call TestCart(AutoCart) |
| 1196 | 1214 |
|
| 1197 | 1215 |
SELECT CASE(NFroz) |
| ... | ... | |
| 2345 | 2363 |
DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local) |
| 2346 | 2364 |
END IF |
| 2347 | 2365 |
|
| 2366 |
If (FPBC) DeAllocate(XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC) |
|
| 2348 | 2367 |
WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt |
| 2349 | 2368 |
|
| 2350 | 2369 |
Close(IOIN) |
| src/Path_module.f90 (revision 10) | ||
|---|---|---|
| 12 | 12 |
INTEGER(KINT), PARAMETER :: MaxFroz=100 |
| 13 | 13 |
REAL(KREAL), PARAMETER :: a0=0.529177249d0 |
| 14 | 14 |
REAL(KREAL), PARAMETER :: Unit=1.d0/a0 |
| 15 |
REAL(KREAL), PARAMETER :: Ang2au=a0, Au2Ang=Unit |
|
| 15 | 16 |
REAL(KREAL) :: Pi |
| 16 | 17 |
|
| 17 | 18 |
! Frozen contains the indices of frozen atoms |
| ... | ... | |
| 192 | 193 |
|
| 193 | 194 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| 194 | 195 |
! |
| 196 |
! Variables for periodic calculations |
|
| 197 |
! |
|
| 198 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 199 |
|
|
| 200 |
! Lattice constants |
|
| 201 |
REAL(KREAL) :: lat_a(3), lat_b(3), lat_c(3) |
|
| 202 |
! Inverse lattice constants |
|
| 203 |
REAL(KREAL) :: Latr(3,3) |
|
| 204 |
! FPBC : True if this is a periodic calculations |
|
| 205 |
LOGICAL :: FPBC |
|
| 206 |
! Possible values for ka, kb, kc in VectorPer (and other periodic operations) |
|
| 207 |
INTEGER(KINT) :: kaBeg,kaEnd,kbBeg,kbEnd,kcBeg,kcEnd |
|
| 208 |
! Number of periodic directions |
|
| 209 |
INTEGER(KINT) :: IPer |
|
| 210 |
! Reference cartesian coordinates |
|
| 211 |
REAL(KREAL), ALLOCATABLE :: XGeomRefPBC(:),YGeomRefPBC(:),ZGeomRefPBC(:) ! Nat |
|
| 212 |
! How shall we print the cartesian coordinates ? |
|
| 213 |
! if V_Direct='DIRECT' then we use fractional coord (ie divided by unit cell |
|
| 214 |
! vectors). |
|
| 215 |
CHARACTER(LCHARS) :: V_direct, V_direct_write |
|
| 216 |
|
|
| 217 |
|
|
| 218 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 219 |
! |
|
| 195 | 220 |
! Variables for VASP input/output |
| 196 | 221 |
! |
| 197 | 222 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| 198 | 223 |
|
| 199 |
REAL(KREAL) :: lat_a(3), lat_b(3), lat_c(3) |
|
| 200 | 224 |
REAL(KREAL), ALLOCATABLE :: X0_vasp(:),Y0_vasp(:), Z0_vasp(:) ! nat |
| 201 |
REAL(KREAL) :: Latr(3,3) |
|
| 202 |
CHARACTER(LCHARS) :: V_direct |
|
| 203 | 225 |
! AutoCart : true if user let PATH determines which atoms should be |
| 204 | 226 |
! described in cartesian when COORD=MIXED |
| 205 | 227 |
LOGICAL :: AutoCart |
| src/NoString.f90 (revision 10) | ||
|---|---|---|
| 1 |
SUBROUTINE noString(String) |
|
| 2 |
|
|
| 3 |
! This subroutine suppress the content of string variable |
|
| 4 |
! strings can be delimited by ' or " |
|
| 5 |
! example: |
|
| 6 |
! Progexe='./siesta', -> Progexe=, |
|
| 7 |
|
|
| 8 |
use VarTypes |
|
| 9 |
|
|
| 10 |
IMPLICIT NONE |
|
| 11 |
|
|
| 12 |
CHARACTER(*), INTENT(INOUT) :: String |
|
| 13 |
CHARACTER(VLCHARS) :: Result |
|
| 14 |
INTEGER(KINT) :: i |
|
| 15 |
CHARACTER(2), PARAMETER :: Del='"' // "'" |
|
| 16 |
|
|
| 17 |
DO WHILE (scan(String,Del)/=0) |
|
| 18 |
I=Scan(String,Del) |
|
| 19 |
Result=String(:I-1) |
|
| 20 |
String=String(I+1:) |
|
| 21 |
I=Scan(String,Del) |
|
| 22 |
String=TRIM(Result) // String(I+1:) |
|
| 23 |
END DO |
|
| 24 |
|
|
| 25 |
END SUBROUTINE NoString |
|
| 26 |
|
|
| src/CalcCnct.f90 (revision 10) | ||
|---|---|---|
| 30 | 30 |
|
| 31 | 31 |
SUBROUTINE CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact) |
| 32 | 32 |
|
| 33 |
use Path_module, only : NMaxL, max_Z,Nom,Prog, KINT, KREAL |
|
| 33 |
use Path_module, only : NMaxL, max_Z,Nom,Prog, KINT, KREAL, FPBC, & |
|
| 34 |
kaBeg,kaEnd,kbBeg,kbEnd,kcBeg,kcEnd |
|
| 34 | 35 |
|
| 35 | 36 |
IMPLICIT NONE |
| 36 | 37 |
|
| ... | ... | |
| 64 | 65 |
END DO |
| 65 | 66 |
|
| 66 | 67 |
if (debug) THEN |
| 67 |
WRITE(*,*) 'CalcCnct' |
|
| 68 |
WRITE(*,*) 'CalcCnct : covalent radii used'
|
|
| 68 | 69 |
DO iat=1,na |
| 69 | 70 |
i=atome(iat) |
| 70 |
WRITE(*,*) Nom(I),I,r_cov(i),r_cov(i)*fact
|
|
| 71 |
WRITE(*,*) Nom(I),I,r_cov(i)*fact |
|
| 71 | 72 |
END DO |
| 73 |
WRITE(*,*) 'Coordinates' |
|
| 74 |
DO iat=1,na |
|
| 75 |
i=atome(iat) |
|
| 76 |
WRITE(*,*) Nom(I),x(iat),y(iat),z(iat) |
|
| 77 |
END DO |
|
| 72 | 78 |
END IF |
| 73 | 79 |
|
| 74 |
IF (PROG/="VASP") THEN
|
|
| 80 |
IF (FPBC) THEN
|
|
| 75 | 81 |
DO i=1,na |
| 76 | 82 |
NbLi=LIAISONS(i,0) |
| 77 | 83 |
DO j=i+1,na |
| 78 |
CALL vecteur(j,i,x,y,z,vx,vy,vz,dist) |
|
| 79 |
dist=dist/fact |
|
| 80 |
distth=(r_cov(atome(i))+r_cov(atome(j)))/100. |
|
| 81 |
! if (debug) WRITE(*,*) atome(i),atome(j),dist,distth |
|
| 82 |
if (dist.le.distth) THEN |
|
| 84 |
Bound=.FALSE. |
|
| 85 |
DO ka=kaBeg,kaEnd |
|
| 86 |
DO Kb=kbBeg,kbEnd |
|
| 87 |
DO Kc=kcBeg,kcEnd |
|
| 88 |
CALL VectorPer(j,i,ka,kb,kc,x,y,z,vx,vy,vz,dist) |
|
| 89 |
dist=dist/fact |
|
| 90 |
distth=(r_cov(atome(i))+r_cov(atome(j)))/100. |
|
| 91 |
! if (debug) WRITE(*,*) atome(i),atome(j),dist,distth |
|
| 92 |
if (dist.le.distth) Bound=.TRUE. |
|
| 93 |
END DO |
|
| 94 |
END DO |
|
| 95 |
END DO |
|
| 96 |
IF (Bound) THEN |
|
| 97 |
if (debug) WRITE(*,*) "Adding a bond between:",i,j |
|
| 83 | 98 |
NbLi=NbLi+1 |
| 84 | 99 |
LIAISONS(i,NbLi)=j; |
| 85 | 100 |
NBlj=LIAISONS(j,0)+1 |
| ... | ... | |
| 93 | 108 |
DO i=1,na |
| 94 | 109 |
NbLi=LIAISONS(i,0) |
| 95 | 110 |
DO j=i+1,na |
| 96 |
Bound=.FALSE. |
|
| 97 |
DO ka=-1,1 |
|
| 98 |
DO Kb=-1,1 |
|
| 99 |
DO Kc=-1,1 |
|
| 100 |
CALL VectorPer(j,i,ka,kb,kc,x,y,z,vx,vy,vz,dist) |
|
| 101 |
dist=dist/fact |
|
| 102 |
distth=(r_cov(atome(i))+r_cov(atome(j)))/100. |
|
| 103 |
! if (debug) WRITE(*,*) atome(i),atome(j),dist,distth |
|
| 104 |
if (dist.le.distth) Bound=.TRUE. |
|
| 105 |
END DO |
|
| 106 |
END DO |
|
| 107 |
END DO |
|
| 108 |
IF (Bound) THEN |
|
| 109 |
if (debug) WRITE(*,*) "Adding a bond between:",i,j |
|
| 111 |
CALL vecteur(j,i,x,y,z,vx,vy,vz,dist) |
|
| 112 |
dist=dist/fact |
|
| 113 |
distth=(r_cov(atome(i))+r_cov(atome(j)))/100. |
|
| 114 |
! if (debug) WRITE(*,*) atome(i),atome(j),dist,distth |
|
| 115 |
if (dist.le.distth) THEN |
|
| 110 | 116 |
NbLi=NbLi+1 |
| 111 | 117 |
LIAISONS(i,NbLi)=j; |
| 112 | 118 |
NBlj=LIAISONS(j,0)+1 |
| src/CheckPeriodicBound.f90 (revision 10) | ||
|---|---|---|
| 1 |
SUBROUTINE CheckPeriodicBound
|
|
| 2 |
! This subroutine is only for periodic systems |
|
| 3 |
! it checks that the atoms do not change their |
|
| 4 |
! positions from one side to the other during the simulations. |
|
| 5 |
! This would ruin the interpolation. |
|
| 1 |
SUBROUTINE CheckPeriodicBound |
|
| 2 |
! This subroutine is only for periodic systems
|
|
| 3 |
! it checks that the atoms do not change their
|
|
| 4 |
! positions from one side to the other during the simulations.
|
|
| 5 |
! This would ruin the interpolation.
|
|
| 6 | 6 |
|
| 7 | 7 |
|
| 8 | 8 |
Use Path_module |
| ... | ... | |
| 15 | 15 |
REAL(KREAL), ALLOCATABLE :: GeomRef(:,:), GeomTmp(:,:) ! (3,Nat) |
| 16 | 16 |
|
| 17 | 17 |
REAL(KREAL) :: Latrloc(3,3),Bloc(3),Xt,Yt,Zt |
| 18 |
REAL(KREAL) :: NormA, NormB,Xrel,Yrel,ZRel,NormC |
|
| 19 |
REAL(KREAL) :: Prod,LatrA(3),LatrB(3) |
|
| 20 |
REAL(KREAL) :: V(3) |
|
| 21 |
|
|
| 18 | 22 |
LOGICAL :: Debug |
| 19 | 23 |
|
| 20 | 24 |
INTERFACE |
| ... | ... | |
| 22 | 26 |
CHARACTER(*), intent(in) :: string |
| 23 | 27 |
logical :: isValid |
| 24 | 28 |
END function VALID |
| 29 |
|
|
| 30 |
|
|
| 31 |
SUBROUTINE die(routine, msg, file, line, unit) |
|
| 32 |
|
|
| 33 |
Use VarTypes |
|
| 34 |
Use io_module |
|
| 35 |
|
|
| 36 |
implicit none |
|
| 37 |
character(len=*), intent(in) :: routine, msg |
|
| 38 |
character(len=*), intent(in), optional :: file |
|
| 39 |
integer(KINT), intent(in), optional :: line, unit |
|
| 40 |
|
|
| 41 |
END SUBROUTINE die |
|
| 42 |
|
|
| 25 | 43 |
END INTERFACE |
| 26 | 44 |
|
| 27 | 45 |
debug=valid('checkperiodicbound')
|
| 28 |
if (debug) WRITE(*,*) "========================= Entering CheckPeriodicBound ==========="
|
|
| 46 |
if (debug) call header("Entering CheckPeriodicBound")
|
|
| 29 | 47 |
|
| 30 | 48 |
ALLOCATE(GeomRef(3,Nat), GeomTmp(3,Nat)) |
| 31 | 49 |
|
| 32 |
V_direct=adjustl(Vasp_direct)
|
|
| 33 |
Call upcase(V_direct)
|
|
| 50 |
! The test is easy in DIRECT space, so we will convert the geometries
|
|
| 51 |
! into this space.
|
|
| 34 | 52 |
|
| 35 |
! The test is easy in DIRECT space, so we will convert the geometries |
|
| 36 |
! into this space. |
|
| 37 |
! We might optimize this as Read_geom converts geometries into NORMAL |
|
| 38 |
! coordinates... so it might be that we are converting them back to direct here... |
|
| 53 |
SELECT CASE (Iper) |
|
| 54 |
CASE(1) |
|
| 55 |
NormA=DOT_PRODUCT(Lat_a,Lat_a) |
|
| 56 |
DO IGeom=2,NGeomI |
|
| 57 |
DO I=1,NAt |
|
| 58 |
DO J=1,3 |
|
| 59 |
v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I) |
|
| 60 |
END DO |
|
| 61 |
XRel=DOT_PRODUCT(Lat_a,v)/NormA |
|
| 62 |
if (abs(Xrel).GE.BoxTol) THEN |
|
| 63 |
If (debug) THEN |
|
| 64 |
WRITE(*,'("Atom ",I5," moved because XRel=",F10.6)') I,XRel
|
|
| 65 |
WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
|
|
| 66 |
WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
|
|
| 67 |
WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Lat_a*Sign(1.0d0,XRel)
|
|
| 68 |
END IF |
|
| 69 |
XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Lat_a*Sign(1.0d0,XRel) |
|
| 70 |
END IF |
|
| 71 |
END DO |
|
| 72 |
END DO |
|
| 73 |
CASE (2) |
|
| 74 |
NormA=DOT_PRODUCT(Lat_a,Lat_a) |
|
| 75 |
NormB=DOT_PRODUCT(Lat_b,Lat_b) |
|
| 76 |
Prod=DOT_PRODUCT(Lat_a,Lat_b) |
|
| 77 |
LatrA=(NormB*Lat_a-Prod*Lat_b)/(NormA*NormB-Prod**2) |
|
| 78 |
LatrB=(NormA*Lat_b-Prod*Lat_a)/(NormA*NormB-Prod**2) |
|
| 39 | 79 |
|
| 80 |
DO IGeom=2,NGeomI |
|
| 81 |
DO I=1,NAt |
|
| 82 |
DO J=1,3 |
|
| 83 |
v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I) |
|
| 84 |
END DO |
|
| 85 |
XRel=DOT_PRODUCT(Latra,v) |
|
| 86 |
YRel=DOT_PRODUCT(Latrb,v) |
|
| 87 |
if (abs(Xrel).GE.BoxTol) THEN |
|
| 88 |
If (debug) THEN |
|
| 89 |
WRITE(*,'("Atom ",I5," moved because XRel=",F10.6)') I,XRel
|
|
| 90 |
WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
|
|
| 91 |
WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
|
|
| 92 |
WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Latra*Sign(1.0d0,XRel)
|
|
| 93 |
END IF |
|
| 94 |
XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Latra*Sign(1.0d0,XRel) |
|
| 95 |
END IF |
|
| 96 |
if (abs(Yrel).GE.BoxTol) THEN |
|
| 97 |
If (debug) THEN |
|
| 98 |
WRITE(*,'("Atom ",I5," moved because YRel=",F10.6)') I,YRel
|
|
| 99 |
WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
|
|
| 100 |
WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
|
|
| 101 |
WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Latrb*Sign(1.0d0,YRel)
|
|
| 102 |
END IF |
|
| 103 |
XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Latrb*Sign(1.0d0,YRel) |
|
| 104 |
END IF |
|
| 105 |
END DO |
|
| 106 |
END DO |
|
| 107 |
|
|
| 108 |
CASE(3) |
|
| 109 |
|
|
| 110 |
|
|
| 111 |
|
|
| 112 |
! if (debug) THEN |
|
| 113 |
! WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_a:",Lat_a |
|
| 114 |
! WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_b:",Lat_b |
|
| 115 |
! WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_c:",Lat_c |
|
| 116 |
! WRITE(*,'(1X,A,3(1X,F10.6))') "BocTol:",BoxTol |
|
| 117 |
|
|
| 118 |
! NormA=DOT_PRODUCT(Lat_a,Lat_a) |
|
| 119 |
! NormB=DOT_PRODUCT(Lat_b,Lat_b) |
|
| 120 |
! Prod=DOT_PRODUCT(Lat_a,Lat_b) |
|
| 121 |
! LatrA=(NormB*Lat_a-Prod*Lat_b)/(NormA*NormB-Prod**2) |
|
| 122 |
! LatrB=(NormA*Lat_b-Prod*Lat_a)/(NormA*NormB-Prod**2) |
|
| 123 |
! ! NormA=DOT_PRODUCT(Latra,Latra) |
|
| 124 |
! ! NormB=DOT_PRODUCT(Latrb,Latrb) |
|
| 125 |
|
|
| 126 |
! DO IGeom=2,NGeomI |
|
| 127 |
! DO I=1,NAt |
|
| 128 |
! DO J=1,3 |
|
| 129 |
! v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I) |
|
| 130 |
! END DO |
|
| 131 |
! XRel=DOT_PRODUCT(Latra,v) |
|
| 132 |
! YRel=DOT_PRODUCT(Latrb,v) |
|
| 133 |
! If (abs(XRel).GT.BoxTol) WRITE(*,*) "Xrel too big:",IGeom,I,Xrel |
|
| 134 |
! If (abs(YRel).GT.BoxTol) WRITE(*,*) "Yrel too big:",IGeom,I,Yrel |
|
| 135 |
! ! WRITE(*,*) "Xrel :",IGeom,I,Xrel |
|
| 136 |
! ! WRITE(*,*) "Yrel :",IGeom,I,Yrel |
|
| 137 |
! ! WRITE(*,*) "Zrel :",IGeom,I,Zrel |
|
| 138 |
|
|
| 139 |
! END DO |
|
| 140 |
! END DO |
|
| 141 |
! END IF |
|
| 142 |
|
|
| 40 | 143 |
Latrloc(1:3,1)=Lat_a |
| 41 | 144 |
Latrloc(1:3,2)=Lat_b |
| 42 | 145 |
Latrloc(1:3,3)=Lat_c |
| 43 | 146 |
Bloc=1. |
| 147 |
! TODO: Replace by LAPACK subroutine |
|
| 44 | 148 |
CALL Gaussj(Latrloc,3,3,Bloc,1,1) |
| 45 | 149 |
|
| 46 |
! We put first geometry as the reference |
|
| 150 |
! We put first geometry as the reference
|
|
| 47 | 151 |
GeomRef=0. |
| 48 | 152 |
DO I=1,Nat |
| 49 | 153 |
DO K=1,3 |
| ... | ... | |
| 53 | 157 |
END DO |
| 54 | 158 |
END DO |
| 55 | 159 |
|
| 56 |
DO IGeom=2,NGeomI |
|
| 57 |
GeomTmp=0. |
|
| 58 |
DO I=1,Nat |
|
| 59 |
DO K=1,3 |
|
| 60 |
DO J=1,3 |
|
| 61 |
GeomTmp(K,I)=GeomTmp(K,I)+XyzGeomI(IGeom,J,I)*LatrLoc(K,J) |
|
| 160 |
DO IGeom=2,NGeomI |
|
| 161 |
GeomTmp=0. |
|
| 162 |
DO I=1,Nat |
|
| 163 |
DO K=1,3 |
|
| 164 |
DO J=1,3 |
|
| 165 |
GeomTmp(K,I)=GeomTmp(K,I)+XyzGeomI(IGeom,J,I)*LatrLoc(K,J) |
|
| 166 |
END DO |
|
| 167 |
if (abs(GeomTmp(K,I)-GeomRef(K,I)).GE.BoxTol) THEN |
|
| 168 |
If (debug) & |
|
| 169 |
WRITE(*,'("Coord ",I1," has been changed for atom ",I5," by ",F4.0," ref:",F9.6," old:",F9.6," new:",F9.6)') &
|
|
| 170 |
K,I,Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)),& |
|
| 171 |
GeomRef(K,I),GeomTmp(K,I), & |
|
| 172 |
GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)) |
|
| 173 |
|
|
| 174 |
GeomTmp(K,I)=GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)) |
|
| 175 |
END IF |
|
| 62 | 176 |
END DO |
| 63 |
if (abs(GeomTmp(K,I)-GeomRef(K,I)).GE.BoxTol) THEN |
|
| 64 |
If (debug) & |
|
| 65 |
WRITE(*,'("Coord ",I1," has been changed for atom ",I5," by ",F4.0," ref:",F9.6," old:",F9.6," new:",F9.6)') &
|
|
| 66 |
K,I,Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)),& |
|
| 67 |
GeomRef(K,I),GeomTmp(K,I), & |
|
| 68 |
GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)) |
|
| 69 |
GeomTmp(K,I)=GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)) |
|
| 70 |
END IF |
|
| 177 |
Xt=GeomTmp(1,I) |
|
| 178 |
Yt=GeomTmp(2,I) |
|
| 179 |
Zt=GeomTmp(3,I) |
|
| 180 |
XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1) |
|
| 181 |
XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2) |
|
| 182 |
XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3) |
|
| 71 | 183 |
END DO |
| 72 |
Xt=GeomTmp(1,I) |
|
| 73 |
Yt=GeomTmp(2,I) |
|
| 74 |
Zt=GeomTmp(3,I) |
|
| 75 |
XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1) |
|
| 76 |
XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2) |
|
| 77 |
XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3) |
|
| 78 | 184 |
END DO |
| 79 |
END DO |
|
| 80 | 185 |
|
| 81 |
if (debug) WRITE(*,*) "========================= Exiting CheckPeriodicBound ===========" |
|
| 186 |
CASE DEFAULT |
|
| 187 |
Call Die("CheckPeriodicBound"," Iper unknown")
|
|
| 188 |
END SELECT |
|
| 189 |
|
|
| 190 |
|
|
| 191 |
if (debug) WRITE(*,*) "========================= Exiting CheckPeriodicBound ===========" |
|
| 82 | 192 |
END SUBROUTINE CheckPeriodicBound |
| 83 | 193 |
|
| 84 | 194 |
|
| src/AlignPartial.f90 (revision 10) | ||
|---|---|---|
| 13 | 13 |
! x2,y2,z2 : coordinates of the second geometry |
| 14 | 14 |
! ListAt : Logical, T if this atom should be considered in the RMSD calc and |
| 15 | 15 |
! alignment. |
| 16 |
! Debug : True if additionnal printing is requiered |
|
| 17 | 16 |
! |
| 18 | 17 |
! output parameters: |
| 19 | 18 |
! U : Real, (3x3) matrix for rotation |
| 20 |
! Rmsd : Read, the rmsd |
|
| 21 | 19 |
! x2,y2,z2 : coordinates of the second geometry ALIGNED to the first one |
| 22 | 20 |
!----------------------------------------------------------------------- |
| 23 | 21 |
|
| src/Egrad.f90 (revision 10) | ||
|---|---|---|
| 8 | 8 |
prog,NCart,XyzGeomF, XyzGeomI, renum, & |
| 9 | 9 |
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & |
| 10 | 10 |
, BprimT,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & |
| 11 |
, OrderInv, XPrimitiveF
|
|
| 11 |
,Order, OrderInv, XPrimitiveF,FPBC,XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC
|
|
| 12 | 12 |
|
| 13 | 13 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) |
| 14 | 14 |
! allocated in Path.f90 |
| ... | ... | |
| 43 | 43 |
REAL(KREAL), ALLOCATABLE :: GradCart(:) |
| 44 | 44 |
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) |
| 45 | 45 |
REAL(KREAL) :: Pi |
| 46 |
REAL(KREAL) :: MRot(3,3), Rmsd |
|
| 46 | 47 |
|
| 47 | 48 |
INTEGER(KINT) :: iat, i, j, IBeg,Idx |
| 48 | 49 |
|
| ... | ... | |
| 92 | 93 |
CASE ('ZMAT')
|
| 93 | 94 |
! In order to avoid problem with numbering and linear angles/diehedral, we convert the |
| 94 | 95 |
! zmat into cartesian coordinates |
| 95 |
! A remplacer par Int2Cart :-) |
|
| 96 | 96 |
Call Int2Cart(nat,indzmat,Geom,GeomCart) |
| 97 |
! In case of PBC calculations, we have to re-orient the geometry onto the user initial geometry |
|
| 98 |
If (FPBC) THEN |
|
| 99 |
! we align this geometry on the initial one |
|
| 100 |
if (debug) THEN |
|
| 101 |
WRITe(*,*) "We are orientating..." |
|
| 102 |
WRITE(*,*) "Geom before orientation" |
|
| 103 |
WRITE(*,*) Nat |
|
| 104 |
WRITE(*,*) "" |
|
| 105 |
DO I=1,Nat |
|
| 106 |
If (renum) Iat=Order(I) |
|
| 107 |
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomCart(I,1),GeomCart(I,2),GeomCart(I,3) |
|
| 108 |
END DO |
|
| 109 |
END IF |
|
| 110 |
|
|
| 111 |
IF (Renum) THEN |
|
| 112 |
DO I=1,Nat |
|
| 113 |
Iat=Order(I) |
|
| 114 |
X(I)=GeomCart(Iat,1) |
|
| 115 |
Y(I)=GeomCart(Iat,2) |
|
| 116 |
Z(I)=GeomCart(Iat,3) |
|
| 117 |
END DO |
|
| 118 |
ELSE |
|
| 119 |
X(1:Nat)=GeomCart(1:Nat,1) |
|
| 120 |
Y(1:Nat)=GeomCart(1:Nat,2) |
|
| 121 |
Z(1:Nat)=GeomCart(1:Nat,3) |
|
| 122 |
END IF |
|
| 123 |
|
|
| 124 |
WRITE(*,*) "Geom before orientation after reorderring" |
|
| 125 |
WRITE(*,*) Nat |
|
| 126 |
WRITE(*,*) "" |
|
| 127 |
DO I=1,Nat |
|
| 128 |
! Iat=I |
|
| 129 |
! If (Renum) Iat=Order(I) |
|
| 130 |
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),X(I),Y(I),Z(I) |
|
| 131 |
END DO |
|
| 132 |
WRITE(*,*) "Ref Geom" |
|
| 133 |
WRITE(*,*) Nat |
|
| 134 |
WRITE(*,*) "" |
|
| 135 |
DO I=1,Nat |
|
| 136 |
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XGeomRefPBC(I),YGeomRefPBC(I),ZGeomRefPBC(I) |
|
| 137 |
END DO |
|
| 138 |
END IF |
|
| 139 |
|
|
| 140 |
Call CalcRmsd(Nat,XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC, & |
|
| 141 |
X,Y,Z, MRot, Rmsd,.TRUE.,.TRUE.) |
|
| 142 |
|
|
| 143 |
If (DEBUG) THEN |
|
| 144 |
WRITE(*,*) "Geom AFTER orientation" |
|
| 145 |
WRITE(*,*) Nat |
|
| 146 |
WRITE(*,*) "" |
|
| 147 |
DO I=1,Nat |
|
| 148 |
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),X(I),Y(I),Z(I) |
|
| 149 |
END DO |
|
| 150 |
END IF |
|
| 151 |
|
|
| 152 |
If (Renum) THEN |
|
| 153 |
Do I=1,Nat |
|
| 154 |
Iat=orderInv(I) |
|
| 155 |
GeomCart(Iat,1)=X(I) |
|
| 156 |
GeomCart(Iat,2)=Y(I) |
|
| 157 |
GeomCart(Iat,3)=Z(I) |
|
| 158 |
END DO |
|
| 159 |
END IF |
|
| 97 | 160 |
CASE ('BAKER')
|
| 98 | 161 |
XPrimRef=XPrimitiveF(IGeom,:) |
| 99 | 162 |
IF (Flag_Opt_Geom) THEN |
| ... | ... | |
| 184 | 247 |
! WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3) |
| 185 | 248 |
!END DO |
| 186 | 249 |
|
| 250 |
|
|
| 251 |
|
|
| 187 | 252 |
SELECT CASE (Prog) |
| 188 | 253 |
CASE ('GAUSSIAN')
|
| 189 | 254 |
! we call the Gaussian routine. |
| src/ReadInput_gaussian.f90 (revision 10) | ||
|---|---|---|
| 13 | 13 |
CHARACTER(*), intent(in) :: string |
| 14 | 14 |
logical :: isValid |
| 15 | 15 |
END function VALID |
| 16 |
|
|
| 17 |
|
|
| 18 |
SUBROUTINE die(routine, msg, file, line, unit) |
|
| 19 |
|
|
| 20 |
Use VarTypes |
|
| 21 |
Use io_module |
|
| 22 |
|
|
| 23 |
implicit none |
|
| 24 |
character(len=*), intent(in) :: routine, msg |
|
| 25 |
character(len=*), intent(in), optional :: file |
|
| 26 |
integer(KINT), intent(in), optional :: line, unit |
|
| 27 |
|
|
| 28 |
END SUBROUTINE die |
|
| 29 |
|
|
| 30 |
|
|
| 16 | 31 |
END INTERFACE |
| 17 | 32 |
|
| 18 | 33 |
|
| 19 |
CHARACTER(132) :: Line |
|
| 34 |
CHARACTER(132) :: Line,LineUp
|
|
| 20 | 35 |
INTEGER(KINT) :: LineL, Idx, Iat |
| 21 |
INTEGER(KINT) :: I |
|
| 36 |
INTEGER(KINT) :: I,NTmp |
|
| 37 |
REAL(KREAL) :: Lat(3,3) |
|
| 22 | 38 |
|
| 23 | 39 |
LOGICAL :: Debug |
| 24 | 40 |
|
| ... | ... | |
| 86 | 102 |
ALLOCATE(Gauss_Paste(NAt)) |
| 87 | 103 |
LineL=1 |
| 88 | 104 |
Iat=0 |
| 105 |
IPer=0 |
|
| 106 |
FPBC=.FALSE. |
|
| 89 | 107 |
Gauss_paste=" " |
| 90 | 108 |
DO While (LineL.GT.0) |
| 91 | 109 |
READ(IOIN,'(A)') Line |
| 92 | 110 |
Line=AdjustL(Line) |
| 93 | 111 |
LineL=Len_TRIM(Line) |
| 94 |
IF (LineL.GT.0) Iat=Iat+1 |
|
| 95 |
Idx=Index(Line,'.',BACK=.TRUE.) |
|
| 96 |
Line=ADJUSTL(Line(Idx+1:)) |
|
| 97 |
|
|
| 98 |
IF (LEN_TRIM(Line).GT.0) THEN |
|
| 99 |
Gauss_paste(Iat)=ADJUSTL(TRIM(Line)) |
|
| 112 |
IF (LineL.GT.0) THEN |
|
| 113 |
LineUp=Line |
|
| 114 |
Call UpCase(LineUp) |
|
| 115 |
if (LineUp(1:2)=="TV") THEN |
|
| 116 |
FPBC=.TRUE. |
|
| 117 |
IPer=IPer+1 |
|
| 118 |
If (Iper>3) THEN |
|
| 119 |
Call Die("ReadInput Gaussian","Iper>3",Unit=IOOUT)
|
|
| 120 |
END IF |
|
| 121 |
NTmp=Index(LineUp," ") |
|
| 122 |
LineUp=LineUp(NTmp:) |
|
| 123 |
Read(LineUp,*) Lat(IPer,1:3) |
|
| 124 |
ELSE |
|
| 125 |
Iat=Iat+1 |
|
| 126 |
! we search for additional information at the end of the line |
|
| 127 |
! for example ONIOM layers |
|
| 128 |
! TO detect the end of the line, we use the fact that all reals |
|
| 129 |
! should contain a '.', and that we have 3 reals/line. |
|
| 130 |
Idx=Index(Line,'.') |
|
| 131 |
Line=ADJUSTL(Line(Idx+1:)) |
|
| 132 |
Idx=Index(Line,'.') |
|
| 133 |
Line=ADJUSTL(Line(Idx+1:)) |
|
| 134 |
Idx=Index(Line,'.') |
|
| 135 |
Line=ADJUSTL(Line(Idx+1:)) |
|
| 136 |
Idx=Index(Line,' ') |
|
| 137 |
If (Idx>0) THEN |
|
| 138 |
Line=ADJUSTL(Line(Idx:)) |
|
| 139 |
IF (LEN_TRIM(Line).GT.0) THEN |
|
| 140 |
Gauss_paste(Iat)=ADJUSTL(TRIM(Line)) |
|
| 141 |
END IF |
|
| 142 |
ELSE |
|
| 143 |
Gauss_paste(Iat)="" |
|
| 144 |
END IF |
|
| 145 |
END IF |
|
| 100 | 146 |
END IF |
| 147 |
|
|
| 101 | 148 |
END DO |
| 102 | 149 |
|
| 150 |
|
|
| 103 | 151 |
IF (Iat.NE.Nat) THEN |
| 104 |
WRITE(*,*) "I found ", Iat," lines for the geometry instead of ",Nat |
|
| 105 |
WRITE(*,*) "ERROR. STOP" |
|
| 106 |
WRITE(IOOUT,*) "I found ", Iat," lines for the geometry instead of ",Nat |
|
| 107 |
WRITE(IOOUT,*) "ERROR. STOP" |
|
| 108 |
STOP |
|
| 152 |
WRITE(Line,*) "I found ", Iat," lines for the geometry instead of ",Nat |
|
| 153 |
Call Die("ReadInput Gaussian","Line",UNIT=IOOUT)
|
|
| 109 | 154 |
END IF |
| 155 |
|
|
| 156 |
IF (FPBC) THEN |
|
| 157 |
Lat_a(1:3)=Lat(1,1:3) |
|
| 158 |
Lat_b(1:3)=Lat(2,1:3) |
|
| 159 |
Lat_c(1:3)=Lat(3,1:3) |
|
| 160 |
If (IPer>=1) THEN |
|
| 161 |
kaBeg=-1 |
|
| 162 |
kaEnd=1 |
|
| 163 |
END IF |
|
| 164 |
If (IPer>=2) THEN |
|
| 165 |
kbBeg=-1 |
|
| 166 |
kbEnd=1 |
|
| 167 |
END IF |
|
| 168 |
If (IPer==3) THEN |
|
| 169 |
kcBeg=-1 |
|
| 170 |
kcEnd=1 |
|
| 171 |
END IF |
|
| 172 |
If (IPer>3) THEN |
|
| 173 |
Call Die("Readinput_gaussian","Found too many Tv lines !",Unit=IOOUT)
|
|
| 174 |
END IF |
|
| 175 |
END IF |
|
| 176 |
|
|
| 177 |
|
|
| 110 | 178 |
! We now read the last part |
| 111 | 179 |
IF (DEBUG) WRITE(*,*) "Reading Gauss End" |
| 112 | 180 |
! READ(IOIN,'(A)') Line |
| ... | ... | |
| 134 | 202 |
END DO |
| 135 | 203 |
|
| 136 | 204 |
WRITE(*,*) |
| 137 |
WRITE(*,*) 'Comment original'
|
|
| 205 |
! WRITE(*,*) '//INFO// Comment original:'
|
|
| 138 | 206 |
|
| 139 | 207 |
Current => Gauss_comment |
| 140 | 208 |
DO WHILE (ASSOCIATED(Current%next)) |
| ... | ... | |
| 146 | 214 |
WRITE(*,*) Trim(Gauss_charge) |
| 147 | 215 |
|
| 148 | 216 |
DO I=1,Nat |
| 149 |
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XyzGeomI(1,1:3,I), TRIM(Gauss_Paste(I)) |
|
| 217 |
WRITE(*,'(1X,A10,3(1X,F15.8),1X,A)') Trim(AtName(I)),XyzGeomI(1,1:3,I), TRIM(Gauss_Paste(I))
|
|
| 150 | 218 |
END DO |
| 151 | 219 |
|
| 152 | 220 |
WRITE(*,*) |
| src/egrad_siesta.f90 (revision 10) | ||
|---|---|---|
| 3 | 3 |
! This routines calculates the energy and the gradient of |
| 4 | 4 |
! a molecule, using Siesta |
| 5 | 5 |
|
| 6 |
use Path_module, only : Nat, renum,Order,OrderInv, Coord, ProgExe |
|
| 6 |
use Path_module, only : Nat, renum,Order,OrderInv, Coord, ProgExe,v_direct,latr
|
|
| 7 | 7 |
use Io_module |
| 8 | 8 |
|
| 9 | 9 |
IMPLICIT NONE |
| ... | ... | |
| 73 | 73 |
logical :: debug, TChk |
| 74 | 74 |
LOGICAL, SAVE :: first=.true. |
| 75 | 75 |
|
| 76 |
REAL(KREAL) :: Pi |
|
| 76 |
REAL(KREAL), ALLOCATABLE :: X(:),Y(:),Z(:) ! Nat |
|
| 77 |
REAL(KREAL), ALLOCATABLE :: GeomCartLoc(:,:) !Nat,3 |
|
| 77 | 78 |
|
| 78 |
INTEGER(KINT) :: iat, i, n3at,Idx
|
|
| 79 |
INTEGER(KINT) :: iat, i, J,Idx
|
|
| 79 | 80 |
|
| 80 | 81 |
CHARACTER(LCHARS) :: FileIn, FileOut, FileForces,FileE |
| 81 | 82 |
|
| ... | ... | |
| 84 | 85 |
|
| 85 | 86 |
! ====================================================================== |
| 86 | 87 |
|
| 87 |
Pi=dacos(-1.0d0) |
|
| 88 |
n3at=3*nat |
|
| 89 |
|
|
| 90 | 88 |
debug=valid('EGRAD').OR.valid('egrad_siesta')
|
| 91 | 89 |
|
| 92 | 90 |
if (debug) Call Header("Entering Egrad_siesta")
|
| ... | ... | |
| 103 | 101 |
END IF |
| 104 | 102 |
END IF |
| 105 | 103 |
|
| 104 |
ALLOCATE(GeomCartLoc(Nat,3)) |
|
| 105 |
|
|
| 106 | 106 |
gradcart=0. |
| 107 | 107 |
E=0. |
| 108 | 108 |
FileIn=Trim(CalcName) // "." // Trim(ISuffix) |
| ... | ... | |
| 122 | 122 |
|
| 123 | 123 |
WRITE(IOTMP,*) |
| 124 | 124 |
|
| 125 |
IF (V_Direct=='DIRECT') THEN |
|
| 126 |
ALLOCATE(X(Nat),Y(Nat),Z(Nat)) |
|
| 127 |
X=0. |
|
| 128 |
Y=0. |
|
| 129 |
Z=0. |
|
| 130 |
DO I=1,Nat |
|
| 131 |
DO J=1,3 |
|
| 132 |
X(I)=X(I)+GeomCart(I,J)*Latr(1,J) |
|
| 133 |
Y(I)=Y(I)+GeomCart(I,J)*Latr(2,J) |
|
| 134 |
Z(I)=Z(I)+GeomCart(I,J)*Latr(3,J) |
|
| 135 |
END DO |
|
| 136 |
END DO |
|
| 137 |
GeomCartLoc(1:Nat,1)=X(1:Nat) |
|
| 138 |
GeomCartLoc(1:Nat,2)=Y(1:Nat) |
|
| 139 |
GeomCartLoc(1:Nat,3)=Z(1:Nat) |
|
| 140 |
DEALLOCATE(X,Y,Z) |
|
| 141 |
ELSE |
|
| 142 |
GeomCartLoc=GeomCart |
|
| 143 |
END IF |
|
| 144 |
|
|
| 125 | 145 |
WRITE(IOTMP,'(1X,A)') '%block AtomicCoordinatesAndAtomicSpecies' |
| 126 | 146 |
|
| 127 | 147 |
DO I=1,Nat |
| ... | ... | |
| 206 | 226 |
WRITE(*,'(1X,3(1X,F15.8))') GradCart(3*I-2:3*I) |
| 207 | 227 |
END DO |
| 208 | 228 |
|
| 229 |
DeAllocate(GeomCartLoc) |
|
| 230 |
|
|
| 209 | 231 |
! Call Die('Egrad_Siesta','Ah ah.',UNIT=IOOUT)
|
| 210 | 232 |
if (debug) Call header("Egrad_siesta Over")
|
| 211 | 233 |
|
| 234 |
|
|
| 235 |
|
|
| 212 | 236 |
RETURN |
| 213 | 237 |
|
| 214 | 238 |
! ====================================================================== |
| src/NoComment.f90 (revision 10) | ||
|---|---|---|
| 1 |
SUBROUTINE noComment(String) |
|
| 2 |
|
|
| 3 |
! This subroutine suppress all comments from a string |
|
| 4 |
! comments can start with ! or # |
|
| 5 |
|
|
| 6 |
IMPLICIT NONE |
|
| 7 |
|
|
| 8 |
integer, parameter :: KINT = kind(1) |
|
| 9 |
|
|
| 10 |
CHARACTER(*), INTENT(INOUT) :: String |
|
| 11 |
INTEGER(KINT) :: i |
|
| 12 |
|
|
| 13 |
I=Scan(String,"!#") |
|
| 14 |
IF (I/=0) THEN |
|
| 15 |
String=String(:I-1) |
|
| 16 |
END IF |
|
| 17 |
|
|
| 18 |
END SUBROUTINE NoComment |
|
| 19 |
|
|
Formats disponibles : Unified diff