Révision 7
| src/ReadInput_siesta.F90 (revision 7) | ||
|---|---|---|
| 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 |
INTERFACE |
|
| 12 |
function valid(string) result (isValid) |
|
| 13 |
CHARACTER(*), intent(in) :: string |
|
| 14 |
logical :: isValid |
|
| 15 |
END function VALID |
|
| 16 |
END INTERFACE |
|
| 17 |
|
|
| 18 |
|
|
| 19 |
CHARACTER(132) :: Line,LineUp,Line2 |
|
| 20 |
INTEGER(KINT) :: LineL, Idx, Iat |
|
| 21 |
INTEGER(KINT) :: IoRead |
|
| 22 |
|
|
| 23 |
LOGICAL :: Debug |
|
| 24 |
LOGICAL :: FSpecies, FCoord |
|
| 25 |
|
|
| 26 |
|
|
| 27 |
Debug=Valid("readinput").OR.Valid("readinput_siesta")
|
|
| 28 |
|
|
| 29 |
if (debug) Call Header("Entering ReadInput_Siesta")
|
|
| 30 |
|
|
| 31 |
! We read the Siesta input file |
|
| 32 |
ALLOCATE(Siesta_Input) |
|
| 33 |
NULLIFY(Siesta_Input%next) |
|
| 34 |
Current => Siesta_Input |
|
| 35 |
|
|
| 36 |
FSpecies=.FALSE. |
|
| 37 |
FCoord=.FALSE. |
|
| 38 |
|
|
| 39 |
READ(IOIN,'(A)',iostat=IoRead) Line |
|
| 40 |
DO WHILE (IoRead==0) |
|
| 41 |
Line=AdjustL(Line) |
|
| 42 |
LineL=len(Trim(Line)) |
|
| 43 |
LineUp=Line |
|
| 44 |
Call UpCase(LineUp) |
|
| 45 |
|
|
| 46 |
! We look for the NumberofAtoms |
|
| 47 |
Idx=INDEX(LineUp,"NUMBEROFATOMS") |
|
| 48 |
IF (Idx/=0) THEN |
|
| 49 |
Idx=Index(Line," ") |
|
| 50 |
Line2=Trim(AdjustL(Line(Idx+1:))) |
|
| 51 |
Read(Line2,*) IAt |
|
| 52 |
if (Iat/=Nat) THEN |
|
| 53 |
Call Die('ReadInput_siesta','Nat in FDF sample different from Nat Path input', &
|
|
| 54 |
__FILE__, __LINE__, IOOUT) |
|
| 55 |
|
|
| 56 |
END IF |
|
| 57 |
END IF |
|
| 58 |
|
|
| 59 |
! We look for the ChemicalSpeciesLabel block |
|
| 60 |
Idx=INDEX(LineUp,"CHEMICALSPECIESLABEL") |
|
| 61 |
IF (Idx/=0) THEN |
|
| 62 |
FSpecies=.TRUE. |
|
| 63 |
END IF |
|
| 64 |
|
|
| 65 |
! We look for SystemLabel |
|
| 66 |
Idx=INDEX(LineUp,"SYSTEMLABEL") |
|
| 67 |
If (Idx/=0) THEN |
|
| 68 |
Idx=Index(Line," ") |
|
| 69 |
Siesta_Label=Trim(adjustl(Line(Idx+1:))) |
|
| 70 |
END IF |
|
| 71 |
|
|
| 72 |
! We look for Coordinates |
|
| 73 |
|
|
| 74 |
IF ((LineL.NE.0).AND.(Idx.EQ.0)) THEN |
|
| 75 |
current%Line=TRIM(Line) |
|
| 76 |
ALLOCATE(current%next) |
|
| 77 |
Current => Current%next |
|
| 78 |
Nullify(Current%next) |
|
| 79 |
END IF |
|
| 80 |
END DO |
|
| 81 |
|
|
| 82 |
if (debug) Call Header("Exiting ReadInput_Siesta")
|
|
| 83 |
|
|
| 84 |
END SUBROUTINE READINPUT_SIESTA |
|
| src/sinangle.f90 (revision 7) | ||
|---|---|---|
| 1 | 1 |
FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2) |
| 2 | 2 |
|
| 3 |
use Path_module, only : Pi,KINT, KREAL
|
|
| 3 |
use Path_module, only : KREAL |
|
| 4 | 4 |
|
| 5 | 5 |
IMPLICIT NONE |
| 6 | 6 |
|
| src/PrintAnaList.f90 (revision 7) | ||
|---|---|---|
| 1 |
SUBROUTINE PrintAnaList(FileUnit) |
|
| 2 |
! This routines prints the list of geometrical variables to monitor |
|
| 3 |
|
|
| 4 |
|
|
| 5 |
use VarTypes |
|
| 6 |
use Path_module |
|
| 7 |
use Io_module |
|
| 8 |
|
|
| 9 |
IMPLICIT NONE |
|
| 10 |
|
|
| 11 |
INTERFACE |
|
| 12 |
function valid(string) result (isValid) |
|
| 13 |
CHARACTER(*), intent(in) :: string |
|
| 14 |
logical :: isValid |
|
| 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 |
|
|
| 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 |
END INTERFACE |
|
| 31 |
! Input |
|
| 32 |
! Unit to print the description |
|
| 33 |
INTEGER(KINT) :: FileUnit |
|
| 34 |
|
|
| 35 |
! Local |
|
| 36 |
LOGICAL :: Debug |
|
| 37 |
INTEGER(KINT) :: I,J,At1,At2,At3,At4 |
|
| 38 |
|
|
| 39 |
Debug=Valid('PrintAnaList')
|
|
| 40 |
|
|
| 41 |
If (Debug) Call Header("Entering PrintAnaList")
|
|
| 42 |
|
|
| 43 |
CurBary => Bary |
|
| 44 |
|
|
| 45 |
DO I=1,NbCom |
|
| 46 |
WRITe(*,*) "Dbg PrintAnalist: I,NbCom",I,NbCom |
|
| 47 |
At1=CurBary%ListAtoms(0) |
|
| 48 |
if (debug) WRITE(*,'("# c ",I4,20(A3,I3))') At1, &
|
|
| 49 |
(AtName(CurBary%ListAtoms(j)),CurBary%ListAtoms(j),j=1,At1) |
|
| 50 |
|
|
| 51 |
WRITE(FileUnit,'("# c ",I4,20(A3,I3))') At1, &
|
|
| 52 |
(AtName(CurBary%ListAtoms(j)),CurBary%ListAtoms(j),j=1,At1) |
|
| 53 |
if (associated(CurBary%Next)) CurBary=> CurBary%Next |
|
| 54 |
END DO |
|
| 55 |
|
|
| 56 |
CurVar => GeomList |
|
| 57 |
DO I=1,NbVar |
|
| 58 |
SELECT CASE (CurVar%Type) |
|
| 59 |
CASE ('BOND')
|
|
| 60 |
! this is a bond |
|
| 61 |
At1=CurVar%At1 |
|
| 62 |
AT2=CurVar%At2 |
|
| 63 |
if (debug) WRITE(*,'("# b ",2(1X,A,I3))') TRIM(AtName(At1)),At1,AtName(At2),At2
|
|
| 64 |
WRITE(FileUnit,'("# b ",2(1X,A,I3))') TRIM(AtName(At1)),At1,AtName(At2),At2
|
|
| 65 |
CASE ('ANGLE')
|
|
| 66 |
! this is an angle |
|
| 67 |
At1=CurVar%At1 |
|
| 68 |
AT2=CurVar%At2 |
|
| 69 |
At3=CurVar%At3 |
|
| 70 |
if (debug) WRITE(*,'("# a ",4(A3,I3))') AtName(At1),At1,AtName(At2),At2,AtName(At3),At3
|
|
| 71 |
WRITE(FileUnit,'("# a ",4(A3,I3))') AtName(At1),At1,AtName(At2),At2,AtName(At3),At3
|
|
| 72 |
CASE ('DIHEDRAL')
|
|
| 73 |
! this is a dihedral |
|
| 74 |
At1=CurVar%At1 |
|
| 75 |
AT2=CurVar%At2 |
|
| 76 |
At3=CurVar%At3 |
|
| 77 |
At4=CurVar%At4 |
|
| 78 |
if (debug) WRITE(*,'("# d ",4(A3,I3))') AtName(At1),At1,AtName(At2),At2,AtName(At3),At3,AtName(At4),At4
|
|
| 79 |
WRITE(FileUnit,'("# d ",4(A3,I3))') AtName(At1),At1,AtName(At2),At2,AtName(At3),At3,AtName(At4),At4
|
|
| 80 |
END SELECT |
|
| 81 |
If (Associated(CurVar%Next)) CurVar => CurVar%next |
|
| 82 |
END DO |
|
| 83 |
|
|
| 84 |
If (Debug) Call Header("Exiting PrintAnaList")
|
|
| 85 |
|
|
| 86 |
END SUBROUTINE PrintAnaList |
|
| src/AnalyzeGeom.f90 (revision 7) | ||
|---|---|---|
| 1 |
SUBROUTINE AnalyzeGeom(GeomCart,Values) |
|
| 2 |
! This routines read a list of geometrical variables to monitor |
|
| 3 |
! This is inspired from Xyz2Path (that was inspired by Xyz2scan ...) |
|
| 4 |
|
|
| 5 |
|
|
| 6 |
use VarTypes |
|
| 7 |
use Path_module |
|
| 8 |
use Io_module |
|
| 9 |
|
|
| 10 |
IMPLICIT NONE |
|
| 11 |
|
|
| 12 |
|
|
| 13 |
INTERFACE |
|
| 14 |
function valid(string) result (isValid) |
|
| 15 |
CHARACTER(*), intent(in) :: string |
|
| 16 |
logical :: isValid |
|
| 17 |
END function VALID |
|
| 18 |
|
|
| 19 |
SUBROUTINE die(routine, msg, file, line, unit) |
|
| 20 |
|
|
| 21 |
Use VarTypes |
|
| 22 |
Use io_module |
|
| 23 |
|
|
| 24 |
implicit none |
|
| 25 |
|
|
| 26 |
character(len=*), intent(in) :: routine, msg |
|
| 27 |
character(len=*), intent(in), optional :: file |
|
| 28 |
integer(KINT), intent(in), optional :: line, unit |
|
| 29 |
|
|
| 30 |
END SUBROUTINE die |
|
| 31 |
|
|
| 32 |
|
|
| 33 |
SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef) |
|
| 34 |
|
|
| 35 |
Use VarTypes |
|
| 36 |
Use Io_module |
|
| 37 |
Use Path_module, only : pi |
|
| 38 |
|
|
| 39 |
IMPLICIT NONE |
|
| 40 |
|
|
| 41 |
Type (ListCoord), POINTER :: Coordinate |
|
| 42 |
INTEGER(KINT), INTENT(IN) :: Nat,NPrim |
|
| 43 |
REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat) |
|
| 44 |
REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) |
|
| 45 |
REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim) |
|
| 46 |
|
|
| 47 |
END SUBROUTINE CALC_XPRIM |
|
| 48 |
|
|
| 49 |
END INTERFACE |
|
| 50 |
|
|
| 51 |
! Input |
|
| 52 |
REAL(KREAL),INTENT(IN) :: GeomCart(Nat,3) |
|
| 53 |
REAL(KREAL), INTENT(OUT) :: Values(NbVar) |
|
| 54 |
|
|
| 55 |
LOGICAL :: Debug |
|
| 56 |
INTEGER(KINT) :: I,J,K,NatT |
|
| 57 |
REAL(KREAL), ALLOCATABLE :: GeoCartLoc(:,:) ! (Nat+NbCom,3) |
|
| 58 |
REAL(KREAL), ALLOCATABLE :: x(:),y(:),z(:) ! Nat+NbCom |
|
| 59 |
REAL(KREAL) :: COG(3),Weight |
|
| 60 |
|
|
| 61 |
|
|
| 62 |
|
|
| 63 |
Debug=Valid('AnaGeom')
|
|
| 64 |
|
|
| 65 |
If (Debug) Call Header("Entering AnalyzeGeom")
|
|
| 66 |
|
|
| 67 |
if (debug) THEN |
|
| 68 |
WRITE(*,*) "AnalyzeGeom - GeomCart" |
|
| 69 |
DO K=1,Nat |
|
| 70 |
WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3) |
|
| 71 |
END DO |
|
| 72 |
END IF |
|
| 73 |
|
|
| 74 |
NAtt=Nat+NbCom |
|
| 75 |
ALLOCATE(GeoCartLoc(Natt,3),x(Natt),y(Natt),z(Natt)) |
|
| 76 |
GeoCartLoc(1:Nat,:)=GeomCart(:,:) |
|
| 77 |
CurBary => Bary |
|
| 78 |
|
|
| 79 |
DO I=1, NbCom |
|
| 80 |
COG=0. |
|
| 81 |
Weight=0. |
|
| 82 |
DO j=1,CurBary%ListAtoms(0) |
|
| 83 |
DO k=1,3 |
|
| 84 |
COG(k)=COG(k)+GeomCart(CurBary%ListAtoms(j),k)*CurBary%Weights(j) |
|
| 85 |
END DO |
|
| 86 |
Weight=Weight+CurBary%Weights(j) |
|
| 87 |
END DO |
|
| 88 |
COG=COG/Weight |
|
| 89 |
DO k=1,3 |
|
| 90 |
GeoCartLoc(Nat+i,k)=COG(k) |
|
| 91 |
END DO |
|
| 92 |
END DO |
|
| 93 |
|
|
| 94 |
Values=0. |
|
| 95 |
|
|
| 96 |
if (debug) THEN |
|
| 97 |
WRITE(*,*) "AnalyzeGeom before Calc_Xprim - GeomCartLoc" |
|
| 98 |
DO K=1,Nat+NbCom |
|
| 99 |
WRITE(*,*) K,GeoCartLoc(K,1:3) |
|
| 100 |
END DO |
|
| 101 |
END IF |
|
| 102 |
|
|
| 103 |
x = GeoCartLoc(:,1) |
|
| 104 |
y = GeoCartLoc(:,2) |
|
| 105 |
z = GeoCartLoc(:,3) |
|
| 106 |
|
|
| 107 |
Call Calc_XPrim(Natt,x,y,z,GeomList,NbVar,Values) |
|
| 108 |
|
|
| 109 |
DeALLOCATE(GeoCartLoc,x,y,z) |
|
| 110 |
|
|
| 111 |
if (debug) THEN |
|
| 112 |
WRITE(*,*) 'AnalyzeGeom: NbVar,Values',NbVar,Values |
|
| 113 |
END IF |
|
| 114 |
|
|
| 115 |
If (Debug) Call Header("Exiting AnalizeGeom")
|
|
| 116 |
|
|
| 117 |
END SUBROUTINE AnalyzeGeom |
|
| src/Io_module.f90 (revision 7) | ||
|---|---|---|
| 8 | 8 |
SAVE |
| 9 | 9 |
|
| 10 | 10 |
INTEGER(KINT) :: IOIN=11, IOOUT=12, IOCART=14 |
| 11 |
INTEGER(KINT) :: IOGEOM=15 |
|
| 11 |
INTEGER(KINT) :: IOGEOM=15, IODAT=16
|
|
| 12 | 12 |
INTEGER(KINT), PARAMETER :: IOTMP=21,IOTMP2=22, IOTMP3=23 |
| 13 | 13 |
INTEGER(KINT), PARAMETER :: IOERR=19 |
| 14 |
|
|
| 15 |
TYPE Input_Line |
|
| 16 |
CHARACTER(LCHARS) :: Line |
|
| 17 |
TYPE (Input_Line), POINTER :: Prev |
|
| 18 |
TYPE (Input_Line), POINTER :: Next |
|
| 19 |
END TYPE Input_Line |
|
| 20 |
|
|
| 21 | 14 |
CHARACTER(SCHARS) :: RunMode |
| 22 | 15 |
|
| 23 | 16 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| src/VarTypes.f90 (revision 7) | ||
|---|---|---|
| 13 | 13 |
INTEGER(KINT) :: At1,At2,At3,At4 |
| 14 | 14 |
! Value of the Coordinate |
| 15 | 15 |
REAL(KREAL) :: Value |
| 16 |
!! ! This is the value to use to multiply the value to print it in a human readable format |
|
| 17 |
!! REAL(KREAL) :: PrintFactor |
|
| 16 | 18 |
INTEGER(KINT) :: SignDihedral |
| 17 | 19 |
TYPE (ListCoord), POINTER :: Next |
| 18 | 20 |
END TYPE ListCoord |
| 19 | 21 |
|
| 22 |
Type Barycenter |
|
| 23 |
INTEGER(KINT), ALLOCATABLE :: ListAtoms(:) ! Nat max |
|
| 24 |
! Weights is not used for now |
|
| 25 |
REAL(KREAL), ALLOCATABLE :: Weights(:) ! Nat max |
|
| 26 |
TYPE (Barycenter), POINTER :: Next |
|
| 27 |
END type BARYCENTER |
|
| 20 | 28 |
|
| 29 |
|
|
| 30 |
TYPE Input_Line |
|
| 31 |
CHARACTER(LCHARS) :: Line |
|
| 32 |
TYPE (Input_Line), POINTER :: Prev |
|
| 33 |
TYPE (Input_Line), POINTER :: Next |
|
| 34 |
END TYPE Input_Line |
|
| 35 |
|
|
| 36 |
|
|
| 21 | 37 |
END MODULE VARTYPES |
| src/Path_module.f90 (revision 7) | ||
|---|---|---|
| 54 | 54 |
! two initial geometries |
| 55 | 55 |
REAL(KREAL) :: BoxTol=0.5 |
| 56 | 56 |
|
| 57 |
|
|
| 58 | 57 |
LOGICAL :: Freq, MW, Bohr, Renum, Hinv |
| 59 | 58 |
! OptReac: if TRUE the reactants will be optimized. Default is FALSE |
| 60 | 59 |
LOGICAL :: OptReac |
| ... | ... | |
| 69 | 68 |
! CalEProd: if TRUE the products energy will be computed. Default is FALSE |
| 70 | 69 |
! Not compatible with RunMode=Para |
| 71 | 70 |
LOGICAL :: CalcEProd |
| 71 |
! AnaGeom: if TRUE the geometries are analyzed |
|
| 72 |
LOGICAL :: AnaGeom |
|
| 73 |
! Name of the file to analyse geometries |
|
| 74 |
CHARACTER(LCHARS) :: AnaFile |
|
| 75 |
! Nb: number of variables to monitor, including Centers of Mass |
|
| 76 |
INTEGER(KINT) :: Nb |
|
| 77 |
! NbVar: number of geometrical variables to monitor |
|
| 78 |
INTEGER(KINT) :: NbVar |
|
| 79 |
! NbCom: number of center of mass to create |
|
| 80 |
INTEGER(KINT) :: NbCom |
|
| 81 |
TYPE(ListCoord), POINTER :: GeomList,CurVar |
|
| 82 |
TYPE(Barycenter), POINTER :: Bary,CurBary |
|
| 83 |
! Format to print the values |
|
| 84 |
CHARACTER(VLCHARS) :: FormAna |
|
| 85 |
! Factor to use to print the values |
|
| 86 |
REAL(KREAL), ALLOCATABLE :: PrintGeomFactor(:) ! NbVar |
|
| 87 |
! How to update the Hessian |
|
| 72 | 88 |
LOGICAL :: IniHup,HupNeighbour |
| 73 | 89 |
|
| 74 | 90 |
REAL(KREAL) :: Fact,FTan |
| src/Extrapol_int.f90 (revision 7) | ||
|---|---|---|
| 15 | 15 |
! Default value of FAlign=.TRUE. |
| 16 | 16 |
use Path_module, only : NMaxPtPath, IntCoordI, Pi, IndZmat, XyzGeomF, & |
| 17 | 17 |
IntCoordF, IntTangent, Renum, Nom, Order, MassAt, SGeom, XyzGeomI, & |
| 18 |
Atome, Nat, NGeomI, NCoord, NGeomF, OrderInv,NFroz,FrozAtoms |
|
| 18 |
Atome, Nat, NGeomI, NCoord, NGeomF, OrderInv,NFroz,FrozAtoms,Align
|
|
| 19 | 19 |
! IndZmat(Nat,5) |
| 20 | 20 |
|
| 21 | 21 |
use Io_module |
| ... | ... | |
| 127 | 127 |
END DO |
| 128 | 128 |
|
| 129 | 129 |
! We align this geometry with the original one |
| 130 |
! PFL 17/July/2006: only if we have more than 4 atoms. |
|
| 131 |
IF (Nat.GE.4) THEN |
|
| 130 |
! PFL 17/July/2006: only if we have more than 4 atoms. |
|
| 131 |
! PFL 2013 Feb 27 ... or if the user asks for it |
|
| 132 |
IF ((Nat.GE.4).OR.Align) THEN |
|
| 132 | 133 |
! The rotation matrix MRot has INTENT(OUT) in subroutine rotation_matrix(...), |
| 133 | 134 |
! which is called in the CalcRmsd(...). |
| 134 | 135 |
! PFL 24 Nov 2008 -> |
| ... | ... | |
| 263 | 264 |
XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3)) |
| 264 | 265 |
END DO |
| 265 | 266 |
|
| 266 |
IF (Nat.GE.4) THEN
|
|
| 267 |
IF ((Nat.GE.4).OR.Align) THEN
|
|
| 267 | 268 |
! PFL 24 Nov 2008 -> |
| 268 | 269 |
! If we have frozen atoms we align only those ones. |
| 269 | 270 |
if (NFroz.GT.0) THEN |
| ... | ... | |
| 288 | 289 |
XYZTmp(I,J)=XyzTMP2(I,J) |
| 289 | 290 |
END DO |
| 290 | 291 |
END DO |
| 292 |
|
|
| 291 | 293 |
|
| 292 | 294 |
|
| 293 | 295 |
s=s+sqrt(ds) |
| ... | ... | |
| 338 | 340 |
WRITE(IOOUT,'(1X,I5)') Nat |
| 339 | 341 |
WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K |
| 340 | 342 |
! PFL 17/July/2006: only if we have more than 4 atoms. |
| 341 |
IF (Nat.GE.4) THEN
|
|
| 343 |
IF ((Nat.GE.4).OR.Align) THEN
|
|
| 342 | 344 |
! PFL 24 Nov 2008 -> |
| 343 | 345 |
! If we have frozen atoms we align only those ones. |
| 344 | 346 |
if (NFroz.GT.0) THEN |
| ... | ... | |
| 382 | 384 |
|
| 383 | 385 |
if (s>=0.9*dist) THEN |
| 384 | 386 |
if (debug) WRITE(*,*) "DBG Extrapol_int L383: Adding last geom" |
| 385 |
write(*,*) "Extrapol_int u,xgeom(NGeomI),s,dist,s-dist",u,xgeom(NGeomI),s,dist,s-dist |
|
| 387 |
if (debug) write(*,*) "Extrapol_int u,xgeom(NGeomI),s,dist,s-dist",u,xgeom(NGeomI),s,dist,s-dist
|
|
| 386 | 388 |
! u=xgeom(NGeomI) |
| 387 | 389 |
s=s-dist |
| 388 | 390 |
|
| ... | ... | |
| 476 | 478 |
WRITE(IOOUT,'(1X,I5)') Nat |
| 477 | 479 |
WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K |
| 478 | 480 |
! PFL 17/July/2006: only if we have more than 4 atoms. |
| 479 |
IF (Nat.GE.4) THEN
|
|
| 481 |
IF ((Nat.GE.4).OR.Align) THEN
|
|
| 480 | 482 |
! PFL 24 Nov 2008 -> |
| 481 | 483 |
! If we have frozen atoms we align only those ones. |
| 482 | 484 |
if (NFroz.GT.0) THEN |
| src/AnaPath.f90 (revision 7) | ||
|---|---|---|
| 1 |
!================================================================ |
|
| 2 |
! |
|
| 3 |
! This subroutine analyzes the geometries along a path |
|
| 4 |
! and prints it to FileUnit |
|
| 5 |
!================================================================ |
|
| 6 |
|
|
| 7 |
SUBROUTINE AnaPath(IOpt,FileUnit) |
|
| 8 |
|
|
| 9 |
|
|
| 10 |
use Io_module |
|
| 11 |
use Path_module, only : Nat, NGeomF, XyzGeomF, AtName,SGeom,MassAt, & |
|
| 12 |
IReparam,NbVar, FormAna, Energies, au2kcal,XyzGeomI,Order,OrderInv, & |
|
| 13 |
Renum, PrintGeomFactor, & |
|
| 14 |
Align, NFroz, FrozAtoms |
|
| 15 |
|
|
| 16 |
|
|
| 17 |
IMPLICIT NONE |
|
| 18 |
|
|
| 19 |
! Input |
|
| 20 |
! Iteration number |
|
| 21 |
INTEGER(KINT), INTENT(IN) :: IOpt |
|
| 22 |
! Unit file to print |
|
| 23 |
INTEGER(KINT), INTENT(IN) :: FileUnit |
|
| 24 |
|
|
| 25 |
|
|
| 26 |
INTERFACE |
|
| 27 |
function valid(string) result (isValid) |
|
| 28 |
CHARACTER(*), intent(in) :: string |
|
| 29 |
logical :: isValid |
|
| 30 |
END function VALID |
|
| 31 |
|
|
| 32 |
FUNCTION test_zmat(na,ind_zmat) |
|
| 33 |
|
|
| 34 |
USE Path_module |
|
| 35 |
|
|
| 36 |
logical :: test_zmat |
|
| 37 |
integer(KINT) :: na |
|
| 38 |
integer(KINT) :: ind_zmat(Na,5) |
|
| 39 |
END FUNCTION TEST_ZMAT |
|
| 40 |
|
|
| 41 |
|
|
| 42 |
SUBROUTINE die(routine, msg, file, line, unit) |
|
| 43 |
|
|
| 44 |
Use VarTypes |
|
| 45 |
Use io_module |
|
| 46 |
|
|
| 47 |
implicit none |
|
| 48 |
|
|
| 49 |
character(len=*), intent(in) :: routine, msg |
|
| 50 |
character(len=*), intent(in), optional :: file |
|
| 51 |
integer(KINT), intent(in), optional :: line, unit |
|
| 52 |
|
|
| 53 |
END SUBROUTINE die |
|
| 54 |
|
|
| 55 |
SUBROUTINE Warning(routine, msg, file, line, unit) |
|
| 56 |
|
|
| 57 |
Use VarTypes |
|
| 58 |
Use io_module |
|
| 59 |
|
|
| 60 |
implicit none |
|
| 61 |
|
|
| 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 Warning |
|
| 67 |
|
|
| 68 |
END INTERFACE |
|
| 69 |
|
|
| 70 |
LOGICAL :: Debug, CalcDist |
|
| 71 |
INTEGER(KINT) :: I,J,K,Iat |
|
| 72 |
REAL(KREAL), ALLOCATABLE :: GeomCart(:,:),GeomCartPrec(:,:) ! Nat,3 |
|
| 73 |
REAL(KREAL), ALLOCATABLE :: xRef(:),yRef(:),zRef(:) ! Nat |
|
| 74 |
REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar |
|
| 75 |
REAL(KREAL) :: ds |
|
| 76 |
REAL(KREAL) :: MRot(3,3), Rmsd |
|
| 77 |
|
|
| 78 |
debug=valid("AnaPath")
|
|
| 79 |
|
|
| 80 |
if (debug) Call header("Entering AnaPath")
|
|
| 81 |
|
|
| 82 |
if (debug) WRITE(*,*) "DBG AnaPaht PrintGeomFactor:",PrintGeomFactor |
|
| 83 |
ALLOCATE(GeomCart(Nat,3)) |
|
| 84 |
ALLOCATE(Values(NbVar)) |
|
| 85 |
|
|
| 86 |
CalcDist=.FALSE. |
|
| 87 |
|
|
| 88 |
WRITE(FileUnit,'("# Iteration ",I5)') Iopt
|
|
| 89 |
|
|
| 90 |
DO J=1,Nat |
|
| 91 |
Iat=J |
|
| 92 |
IF (renum) Iat=Order(J) |
|
| 93 |
GeomCart(J,:)=XyzGeomF(1,:,Iat) |
|
| 94 |
END DO |
|
| 95 |
|
|
| 96 |
|
|
| 97 |
! Do we know the curvilinear values ? |
|
| 98 |
IF (MOD(IOpt,IReparam)/=0) THEN |
|
| 99 |
CalcDist=.TRUE. |
|
| 100 |
SGeom=0. |
|
| 101 |
ALLOCATE(xRef(Nat),yRef(Nat),zref(Nat),GeomCartPrec(Nat,3)) |
|
| 102 |
xRef=GeomCart(:,1) |
|
| 103 |
yRef=GeomCart(:,2) |
|
| 104 |
zRef=GeomCart(:,3) |
|
| 105 |
GeomCartPrec=GeomCart |
|
| 106 |
END IF |
|
| 107 |
|
|
| 108 |
if (debug) THEN |
|
| 109 |
WRITE(*,*) "AnaPath - GeomCart - I=",1 |
|
| 110 |
DO K=1,Nat |
|
| 111 |
WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3) |
|
| 112 |
END DO |
|
| 113 |
WRITE(*,*) "AnaPath - XYzGeomF(I,:,:) - I=",1 |
|
| 114 |
DO K=1,Nat |
|
| 115 |
WRITE(*,'(1X,I5,3(1X,F15.8))') K,XyzGeomF(1,:,K) |
|
| 116 |
END DO |
|
| 117 |
END IF |
|
| 118 |
|
|
| 119 |
Call AnalyzeGeom(GeomCart,Values) |
|
| 120 |
WRITE(FileUnit,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1) |
|
| 121 |
if (debug) WRITE(*,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1) |
|
| 122 |
|
|
| 123 |
DO I=2,NGeomF |
|
| 124 |
DO J=1,Nat |
|
| 125 |
Iat=J |
|
| 126 |
IF (renum) Iat=Order(J) |
|
| 127 |
GeomCart(J,:)=XyzGeomF(I,:,Iat) |
|
| 128 |
END DO |
|
| 129 |
|
|
| 130 |
If (CalcDist) THEN |
|
| 131 |
! PFL 2013 Feb 26 |
|
| 132 |
! For now, I do _NOT_ align the geometries |
|
| 133 |
if (debug) THEN |
|
| 134 |
WRITE(*,*) "AnaPath - GeomCart avant align - I=",I |
|
| 135 |
DO K=1,Nat |
|
| 136 |
WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3) |
|
| 137 |
END DO |
|
| 138 |
END IF |
|
| 139 |
|
|
| 140 |
IF ((Nat.GE.4).OR.Align) THEN |
|
| 141 |
WRITE(*,*) "DBG AnaPath: Aglin 2 geoms" |
|
| 142 |
|
|
| 143 |
! If we have frozen atoms we align only those ones. |
|
| 144 |
if (NFroz.GT.0) THEN |
|
| 145 |
Call AlignPartial(Nat,xRef,yRef,zRef, & |
|
| 146 |
GeomCart(1,1),GeomCart(1,2),GeomCart(1,3), & |
|
| 147 |
FrozAtoms,MRot) |
|
| 148 |
ELSE |
|
| 149 |
Call CalcRmsd(Nat, xRef,yRef,zRef, & |
|
| 150 |
GeomCart(1,1),GeomCart(1,2),GeomCart(1,3), & |
|
| 151 |
MRot,rmsd,.TRUE.,.TRUE.) |
|
| 152 |
END IF |
|
| 153 |
! <- PFL 24 Nov 2008 |
|
| 154 |
|
|
| 155 |
END IF |
|
| 156 |
|
|
| 157 |
if (debug) WRITE(*,*) "Mass:",MassAt(1:Nat) |
|
| 158 |
ds=0. |
|
| 159 |
DO J=1,Nat |
|
| 160 |
Iat=J |
|
| 161 |
! IF (renum) Iat=Order(J) |
|
| 162 |
if (debug) WRITE(*,*) "Mass(J):",J,MassAt(Iat) |
|
| 163 |
|
|
| 164 |
do K=1,3 |
|
| 165 |
ds=ds+MassAt(Iat)*(GeomCartPrec(J,K)-GeomCart(J,K))**2 |
|
| 166 |
END DO |
|
| 167 |
if (debug) WRITE(*,*) "DBG DS:",J,ds |
|
| 168 |
END DO |
|
| 169 |
ds=sqrt(ds) |
|
| 170 |
if (debug) THEN |
|
| 171 |
WRITE(*,*) "AnaPath - GeomCartPrec - I=",I |
|
| 172 |
DO K=1,Nat |
|
| 173 |
WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCartPrec(K,1:3) |
|
| 174 |
END DO |
|
| 175 |
WRITE(*,*) "DBG Anapaht, ds=",ds |
|
| 176 |
END IF |
|
| 177 |
|
|
| 178 |
SGeom(I)=SGeom(I-1)+ds |
|
| 179 |
GeomCartPrec=GeomCart |
|
| 180 |
xRef=GeomCart(:,1) |
|
| 181 |
yRef=GeomCart(:,2) |
|
| 182 |
zRef=GeomCart(:,3) |
|
| 183 |
|
|
| 184 |
END IF |
|
| 185 |
|
|
| 186 |
if (debug) THEN |
|
| 187 |
WRITE(*,*) "AnaPath - GeomCart - I=",I |
|
| 188 |
DO K=1,Nat |
|
| 189 |
WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3) |
|
| 190 |
END DO |
|
| 191 |
END IF |
|
| 192 |
|
|
| 193 |
Call AnalyzeGeom(GeomCart,Values) |
|
| 194 |
WRITE(FileUnit,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*au2kcal,Energies(i) |
|
| 195 |
END DO |
|
| 196 |
|
|
| 197 |
WRITE(FileUnit,*) "" |
|
| 198 |
WRITE(FileUnit,*) "" |
|
| 199 |
|
|
| 200 |
DeAllocate(GeomCart,Values) |
|
| 201 |
If (CalcDist) Deallocate(GeomCartPrec,xRef,yRef,zRef) |
|
| 202 |
|
|
| 203 |
if (debug) Call header("AnaPath over")
|
|
| 204 |
|
|
| 205 |
END SUBROUTINE AnaPath |
|
| 206 |
|
|
| src/Makefile (revision 7) | ||
|---|---|---|
| 232 | 232 |
ReadInput_vasp.f90 \ |
| 233 | 233 |
ReadInput_mopac.f90 \ |
| 234 | 234 |
ReadInput_siesta.f90 \ |
| 235 |
ReadAnaList.f90 \ |
|
| 236 |
PrintAnaList.f90 \ |
|
| 237 |
AnaPath.f90 \ |
|
| 238 |
AnalyzeGeom.f90 \ |
|
| 235 | 239 |
Calc_zmat.f90 \ |
| 236 | 240 |
Calc_zmat_frag.f90 \ |
| 237 | 241 |
Calc_zmat_constr_frag.f90 \ |
| src/Calc_Xprim.f90 (revision 7) | ||
|---|---|---|
| 80 | 80 |
|
| 81 | 81 |
debug=valid("Calc_Xprim")
|
| 82 | 82 |
debugPFL=valid("BakerPFL")
|
| 83 |
if (debug) WRITE(*,*) "============= Entering Cal_XPrim ==============" |
|
| 84 | 83 |
|
| 84 |
if (debug) Call Header("Entering Cal_XPrim")
|
|
| 85 | 85 |
|
| 86 |
|
|
| 86 | 87 |
IF (debug) THEN |
| 87 | 88 |
WRITE(*,*) "DBG Calc_Xprim, geom" |
| 88 | 89 |
DO I=1,Nat |
| 89 | 90 |
WRITE(*,'(1X,I5,3(1X,F15.3))') I,X(I),y(i),z(i) |
| 90 | 91 |
END DO |
| 91 |
WRITE(*,*) "XPrimRef" |
|
| 92 |
WRITE(*,'(15(1X,F10.6))') XPrimRef |
|
| 92 |
IF (Present(XPrimRef)) THEN |
|
| 93 |
WRITE(*,*) "XPrimRef" |
|
| 94 |
WRITE(*,'(15(1X,F10.6))') XPrimRef |
|
| 95 |
END IF |
|
| 96 |
WRITE(*,*) "NPrim:",NPrim |
|
| 93 | 97 |
END IF |
| 94 | 98 |
|
| 95 | 99 |
|
| 96 |
! WRITE(*,*) "Coordinate:",associated(Coordinate)
|
|
| 100 |
WRITE(*,*) "Coordinate:",associated(Coordinate),coordinate%Type
|
|
| 97 | 101 |
|
| 98 |
ScanCoord=>Coordinate
|
|
| 102 |
ScanCoord => Coordinate
|
|
| 99 | 103 |
|
| 100 |
! WRITE(*,*) "ScanCoord:",associated(ScanCoord),ScanCoord%Type
|
|
| 104 |
WRITE(*,*) "ScanCoord:",associated(ScanCoord),ScanCoord%Type |
|
| 101 | 105 |
|
| 102 | 106 |
! WRITE(*,*) "coucou" |
| 103 | 107 |
I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). |
| ... | ... | |
| 120 | 124 |
vx4,vy4,vz4,norm4) |
| 121 | 125 |
Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, & |
| 122 | 126 |
vx5,vy5,vz5,norm5) |
| 123 |
|
|
| 124 |
DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
|
|
| 125 |
vx2,vy2,vz2,norm2)
|
|
| 126 |
Xprimitive(I) = DiheTmp*Pi/180.
|
|
| 127 |
|
|
| 128 |
DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
|
| 129 |
vx2,vy2,vz2,norm2) |
|
| 130 |
Xprimitive(I) = DiheTmp*Pi/180. |
|
| 127 | 131 |
! PFL 15th March 2008 -> |
| 128 | 132 |
! I think that the test on changes less than Pi should be enough |
| 129 | 133 |
!! We treat large dihedral angles differently as +180=-180 mathematically and physically |
| ... | ... | |
| 139 | 143 |
!!!! <- PFL 15th March 2008 |
| 140 | 144 |
! Another test... will all this be consistent ??? |
| 141 | 145 |
! We want the shortest path, so we check that the change in dihedral angles is less than Pi: |
| 142 |
IF (Present(XPrimRef)) THEN |
|
| 143 |
IF (abs(Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN |
|
| 144 |
if (debug) THEN |
|
| 145 |
WRITE(*,*) "Pb dihedral, i,Xprimivite,XPrimref=",i,XPrimitive(I),XPrimRef(I) |
|
| 146 |
WRITE(*,*) "In deg Xprimivite,XPrimref=",XPrimitive(I)*180./Pi,XPrimRef(I)*180/Pi |
|
| 147 |
END IF |
|
| 148 |
if ((Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN |
|
| 149 |
Xprimitive(I)=Xprimitive(I)-2.d0*Pi |
|
| 150 |
ELSE |
|
| 151 |
Xprimitive(I)=Xprimitive(I)+2.d0*Pi |
|
| 152 |
END IF |
|
| 153 |
END IF |
|
| 154 |
if (debug) WRITE(*,*) " New Xprimivite",XPrimitive(I),XPrimitive(I)*180./Pi |
|
| 155 |
END IF |
|
| 156 |
END SELECT |
|
| 146 |
IF (Present(XPrimRef)) THEN |
|
| 147 |
IF (abs(Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN |
|
| 148 |
if (debug) THEN |
|
| 149 |
WRITE(*,*) "Pb dihedral, i,Xprimivite,XPrimref=",i,XPrimitive(I),XPrimRef(I) |
|
| 150 |
WRITE(*,*) "In deg Xprimivite,XPrimref=",XPrimitive(I)*180./Pi,XPrimRef(I)*180/Pi |
|
| 151 |
END IF |
|
| 152 |
if ((Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN |
|
| 153 |
Xprimitive(I)=Xprimitive(I)-2.d0*Pi |
|
| 154 |
ELSE |
|
| 155 |
Xprimitive(I)=Xprimitive(I)+2.d0*Pi |
|
| 156 |
END IF |
|
| 157 |
END IF |
|
| 158 |
if (debug) WRITE(*,*) " New Xprimivite",XPrimitive(I),XPrimitive(I)*180./Pi |
|
| 159 |
END IF |
|
| 160 |
CASE DEFAULT |
|
| 161 |
WRITE(*,*) "Scancoord%type unknown for I",I,scancoord%type |
|
| 162 |
END SELECT |
|
| 157 | 163 |
ScanCoord => ScanCoord%next |
| 158 | 164 |
END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
| 159 | 165 |
|
| ... | ... | |
| 161 | 167 |
|
| 162 | 168 |
IF (debug) THEN |
| 163 | 169 |
WRITE(*,*) "DBG Calc_Xprim Values" |
| 170 |
WRITE(*,*) "Found ",I," primitives" |
|
| 164 | 171 |
|
| 165 | 172 |
ScanCoord=>Coordinate |
| 166 | 173 |
I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). |
| ... | ... | |
| 179 | 186 |
ScanCoord => ScanCoord%next |
| 180 | 187 |
END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
| 181 | 188 |
END IF |
| 182 |
if (debug) WRITE(*,*) "============= Cal_XPrim OVER ==============" |
|
| 189 |
|
|
| 190 |
if (debug) Call Header(" Cal_XPrim OVER ")
|
|
| 191 |
|
|
| 183 | 192 |
END SUBROUTINE Calc_Xprim |
| src/bib_oper.f90 (revision 7) | ||
|---|---|---|
| 104 | 104 |
v2x,v2y,v2z, & |
| 105 | 105 |
v3x,v3y,v3z,norm3) |
| 106 | 106 |
|
| 107 |
use Path_module, only : Pi,KINT, KREAL
|
|
| 107 |
use Path_module, only : KREAL |
|
| 108 | 108 |
|
| 109 | 109 |
real(KREAL) :: v1x,v1y,v1z ! what do you do with norm1, norm2??? |
| 110 | 110 |
real(KREAL) :: v2x,v2y,v2z |
| src/WriteInput_siesta.f90 (revision 7) | ||
|---|---|---|
| 1 |
SUBROUTINE ReadInput_siesta
|
|
| 1 |
SUBROUTINE WriteInput_siesta(GeomCart,FileUnit)
|
|
| 2 | 2 |
|
| 3 |
! This routine reads an input template for Siesta
|
|
| 3 |
! This routine writes an input for Siesta
|
|
| 4 | 4 |
|
| 5 | 5 |
use VarTypes |
| 6 | 6 |
use Path_module |
| ... | ... | |
| 14 | 14 |
logical :: isValid |
| 15 | 15 |
END function VALID |
| 16 | 16 |
|
| 17 |
FUNCTION SearchInput(Input,String,Line,Clean) Result (Found) |
|
| 18 | 17 |
|
| 19 |
Use Vartypes |
|
| 20 |
use io_module |
|
| 21 |
! Input |
|
| 22 |
TYPE (Input_line), POINTER, INTENT(IN) :: Input |
|
| 23 |
CHARACTER(*), INTENT(IN) :: String |
|
| 24 |
CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean |
|
| 18 |
SUBROUTINE WriteList(Input,Unit) |
|
| 19 |
! This routine reads an input template for Siesta |
|
| 25 | 20 |
|
| 26 |
! Output
|
|
| 27 |
TYPE (Input_line), POINTER, INTENT(OUT) :: Line
|
|
| 21 |
use VarTypes
|
|
| 22 |
use Io_module
|
|
| 28 | 23 |
|
| 29 |
LOGICAL :: Found
|
|
| 24 |
IMPLICIT NONE
|
|
| 30 | 25 |
|
| 31 |
END FUNCTION SearchInput |
|
| 26 |
! Input |
|
| 27 |
TYPE (Input_line), POINTER, INTENT(IN) :: Input |
|
| 28 |
INTEGER(KINT), OPTIONAL, INTENT(IN) :: Unit |
|
| 29 |
|
|
| 30 |
END SUBROUTINE WriteList |
|
| 32 | 31 |
|
| 33 |
SUBROUTINE die(routine, msg, file, line, unit) |
|
| 34 |
|
|
| 35 |
Use VarTypes |
|
| 36 |
Use io_module |
|
| 37 |
|
|
| 38 |
implicit none |
|
| 39 |
!--------------------------------------------------------------- Input Variables |
|
| 40 |
character(len=*), intent(in) :: routine, msg |
|
| 41 |
character(len=*), intent(in), optional :: file |
|
| 42 |
integer(KINT), intent(in), optional :: line, unit |
|
| 43 |
|
|
| 44 |
END SUBROUTINE die |
|
| 45 |
|
|
| 46 |
SUBROUTINE Warning(routine, msg, file, line, unit) |
|
| 47 |
|
|
| 48 |
Use VarTypes |
|
| 49 |
Use io_module |
|
| 50 |
|
|
| 51 |
implicit none |
|
| 52 |
!--------------------------------------------------------------- Input Variables |
|
| 53 |
character(len=*), intent(in) :: routine, msg |
|
| 54 |
character(len=*), intent(in), optional :: file |
|
| 55 |
integer(KINT), intent(in), optional :: line, unit |
|
| 56 |
|
|
| 57 |
END SUBROUTINE Warning |
|
| 58 |
|
|
| 59 | 32 |
END INTERFACE |
| 60 | 33 |
|
| 61 |
|
|
| 62 |
CHARACTER(132) :: Line,LineUp,Line2 |
|
| 63 |
INTEGER(KINT) :: LineL, Idx |
|
| 64 |
INTEGER(KINT) :: ISpec, IAt,I |
|
| 65 |
INTEGER(KINT) :: IoRead, ITmp, JTmp |
|
| 66 |
REAL(KREAL) :: Xtmp, Ytmp, Ztmp |
|
| 34 |
!Input |
|
| 35 |
! Geometry in cartesian coordinates |
|
| 36 |
REAL(KREAL), INTENT (IN) :: geomcart(Nat,3) |
|
| 37 |
! Unit to write to |
|
| 38 |
INTEGER(KINT), INTENT(IN) :: FileUnit |
|
| 67 | 39 |
|
| 68 | 40 |
LOGICAL :: Debug |
| 69 |
LOGICAL :: FSpecies, FCoord
|
|
| 41 |
INTEGER(KINT) :: I,Iat
|
|
| 70 | 42 |
|
| 71 |
TYPE(Input_Line), POINTER :: Search,Bla
|
|
| 43 |
Debug=Valid("WriteInput").OR.Valid("WriteInput_siesta")
|
|
| 72 | 44 |
|
| 73 |
Debug=Valid("readinput").OR.Valid("readinput_siesta")
|
|
| 45 |
if (debug) Call Header("Entering WriteInput_Siesta")
|
|
| 74 | 46 |
|
| 75 |
if (debug) Call Header("Entering ReadInput_Siesta")
|
|
| 76 | 47 |
|
| 77 |
! We read the Siesta input file
|
|
| 48 |
Call WriteList(Siesta_input,Unit=FileUnit)
|
|
| 78 | 49 |
|
| 79 |
ALLOCATE(Siesta_Input) |
|
| 80 |
NULLIFY(Siesta_Input%next) |
|
| 81 |
NULLIFY(Siesta_Input%prev) |
|
| 82 |
Current => Siesta_Input |
|
| 50 |
WRITE(FileUnit,*) |
|
| 83 | 51 |
|
| 84 |
FSpecies=.FALSE. |
|
| 85 |
FCoord=.FALSE. |
|
| 52 |
WRITE(FileUnit,'(1X,A)') '%block AtomicCoordinatesAndAtomicSpecies' |
|
| 86 | 53 |
|
| 87 |
READ(IOIN,'(A)',iostat=IoRead) Line |
|
| 54 |
DO I=1,Nat |
|
| 55 |
If (renum) THEN |
|
| 56 |
Iat=Order(I) |
|
| 57 |
WRITE(FileUnit,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(Iat,:)/Siesta_Unit_Read, IdxSpecies(Iat),TRIM(Siesta_Paste(I)) |
|
| 58 |
ELSE |
|
| 59 |
Iat=OrderInv(I) |
|
| 60 |
WRITE(FileUnit,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(I,:)/Siesta_Unit_Read, IdxSpecies(Iat), TRIM(Siesta_Paste(Iat)) |
|
| 61 |
END IF |
|
| 62 |
END DO |
|
| 88 | 63 |
|
| 89 |
DO WHILE (IoRead==0) |
|
| 90 |
Line=AdjustL(Line) |
|
| 91 |
if (debug) WRITE(*,*) 'Line:', Line |
|
| 92 |
current%Line=TRIM(Line) |
|
| 93 |
ALLOCATE(current%next) |
|
| 94 |
Current%next%prev => Current |
|
| 95 |
Current => Current%next |
|
| 96 |
Nullify(Current%next) |
|
| 64 |
WRITE(FileUnit,'(1X,A)') '%endblock AtomicCoordinatesAndAtomicSpecies' |
|
| 65 |
WRITE(FileUnit,*) |
|
| 97 | 66 |
|
| 98 |
READ(IOIN,'(A)',iostat=IoRead) Line |
|
| 99 |
END DO |
|
| 100 | 67 |
|
| 101 |
! We analyse the input |
|
| 102 | 68 |
|
| 103 |
! We look for the NumberofAtoms |
|
| 104 |
If (SearchInput(Siesta_Input,"NUMBEROFATOMS",Search,Clean=".-_")) THEN |
|
| 105 |
Line=AdjustL(Search%line) |
|
| 106 |
Idx=Index(Line," ") |
|
| 107 |
Line2=Trim(AdjustL(Line(Idx+1:))) |
|
| 108 |
Read(Line2,*) IAt |
|
| 109 |
if (Iat/=Nat) THEN |
|
| 110 |
Call Die('ReadInput_siesta','Nat in FDF sample different from Nat Path input', Unit=IOOUT)
|
|
| 111 |
END IF |
|
| 112 |
ELSE |
|
| 113 |
! There is no atom number defined !!! |
|
| 114 |
Call Warning('ReadInput_siesta','No NumberOfAtoms in FDF sample', Unit=IOOUT)
|
|
| 115 |
WRITE(current%Line,'(1X,"NumberOfAtoms ",I5)') Nat |
|
| 116 |
ALLOCATE(current%next) |
|
| 117 |
Current%next%prev => Current |
|
| 118 |
Current => Current%next |
|
| 119 |
Nullify(Current%next) |
|
| 120 |
|
|
| 121 |
END IF |
|
| 122 |
|
|
| 123 |
! We look for the NumberofSpecies |
|
| 124 |
IF (SearchInput(Siesta_Input,"NUMBEROFSPECIES",Search,Clean=".-_")) THEN |
|
| 125 |
Line=AdjustL(Search%line) |
|
| 126 |
Idx=Index(Line," ") |
|
| 127 |
Line2=Trim(AdjustL(Line(Idx+1:))) |
|
| 128 |
Read(Line2,*) Siesta_NbSpecies |
|
| 129 |
END IF |
|
| 130 |
|
|
| 131 |
! We look for SystemLabel |
|
| 132 |
If (SearchInput(Siesta_Input,"SYSTEMLABEL",Search,Clean=".-_")) THEN |
|
| 133 |
Line=AdjustL(Search%line) |
|
| 134 |
Idx=Index(Line," ") |
|
| 135 |
Siesta_Label=Trim(adjustl(Line(Idx+1:))) |
|
| 136 |
ELSE |
|
| 137 |
Siesta_label='siesta' |
|
| 138 |
END IF |
|
| 139 |
|
|
| 140 |
! We look for the ChemicalSpeciesLabel block |
|
| 141 |
IF (SearchInput(Siesta_Input,"CHEMICALSPECIESLABEL",Search,Clean=".-_")) THEN |
|
| 142 |
ALLOCATE(ListSpecies(Siesta_NbSpecies)) |
|
| 143 |
ALLOCATE(Siesta_SpeciesName(Siesta_NbSpecies)) |
|
| 144 |
ALLOCATE(IdxSpecies(NAt)) |
|
| 145 |
DO I=1,Siesta_NbSpecies |
|
| 146 |
Search => Search%next |
|
| 147 |
Line=AdjustL(Search%Line) |
|
| 148 |
Read(Line,*) ITmp |
|
| 149 |
Idx=Index(Line,' ') |
|
| 150 |
Line=AdjustL(Line(Idx+1:)) |
|
| 151 |
Read(Line,*) Ztmp |
|
| 152 |
ListSpecies(ITmp)=ZTmp |
|
| 153 |
Idx=Index(Line,' ') |
|
| 154 |
Line=AdjustL(Line(Idx+1:)) |
|
| 155 |
Idx=Index(Line,' ') |
|
| 156 |
Siesta_SpeciesName(ITmp)=Line(1:Idx-1) |
|
| 157 |
END DO |
|
| 158 |
if (Debug) THEN |
|
| 159 |
WRITE(*,*) 'Found ',Siesta_NbSpecies,' species' |
|
| 160 |
DO I=1,Siesta_NbSpecies |
|
| 161 |
WRITE(*,*) I, ListSpecies(I),TRIM(Siesta_speciesName(I)) |
|
| 162 |
END DO |
|
| 163 |
END IF |
|
| 164 |
|
|
| 165 |
|
|
| 166 |
! We look for the AtomicCoordinatesAndAtomicSpecies block |
|
| 167 |
IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESANDATOMICSPECIES",Search,Clean=".-_")) THEN |
|
| 168 |
ALLOCATE(Siesta_Paste(Nat)) |
|
| 169 |
Current=>Search |
|
| 170 |
DO I=1,NAt |
|
| 171 |
Search => Search%next |
|
| 172 |
Read(Search%line,*) XTmp,YTmp,ZTmp,IdxSpecies(I) |
|
| 173 |
! We save everything but the x,y,z, and species description |
|
| 174 |
Line=AdjustL(search%line) |
|
| 175 |
! We skip x |
|
| 176 |
Idx=Index(Line,' ') |
|
| 177 |
Line=AdjustL(Line(Idx+1:)) |
|
| 178 |
! we skip y |
|
| 179 |
Idx=Index(Line,' ') |
|
| 180 |
Line=AdjustL(Line(Idx+1:)) |
|
| 181 |
! we skip z |
|
| 182 |
Idx=Index(Line,' ') |
|
| 183 |
Line=AdjustL(Line(Idx+1:)) |
|
| 184 |
! we skip species |
|
| 185 |
Idx=Index(Line,' ') |
|
| 186 |
Line=AdjustL(Line(Idx+1:)) |
|
| 187 |
Siesta_Paste(I)=TRIM(Line) |
|
| 188 |
END DO |
|
| 189 |
! We will now delete this block from our input sample as it will then be |
|
| 190 |
! written directly by Opt'n Path |
|
| 191 |
! Search%next point on %endblock |
|
| 192 |
search => search%next |
|
| 193 |
IF (ASSOCIATED(Search%next)) THEN |
|
| 194 |
! we are not at the end of the input file |
|
| 195 |
if (ASSOCIATED(Current%Prev)) THEN |
|
| 196 |
Search%next%prev => current%prev |
|
| 197 |
ELSE |
|
| 198 |
! the coordinate block is at the begining of the input |
|
| 199 |
siesta_input => Search%next |
|
| 200 |
Nullify(Siesta_Input%Prev) |
|
| 201 |
END IF |
|
| 202 |
ELSE |
|
| 203 |
! the coordinate block is the last part of the input |
|
| 204 |
nullify(current%prev%next) |
|
| 205 |
END IF |
|
| 206 |
DO I=1,Nat+2 |
|
| 207 |
bla=>current |
|
| 208 |
current=> current%next |
|
| 209 |
deallocate(bla) |
|
| 210 |
END DO |
|
| 211 |
ELSE |
|
| 212 |
IF (SearchInput(Siesta_Input,"ZMATRIX",Search,Clean=".-_")) THEN |
|
| 213 |
call Die('ReadInput_Siesta','For now, I need the full block'//&
|
|
| 214 |
'AtomicCoordinatesAndAtomicSpecies, ZMatrix not yet handled.') |
|
| 215 |
ELSE |
|
| 216 |
call Die('ReadInput_Siesta','For now, I need the full block' //&
|
|
| 217 |
'AtomicCoordinatesAndAtomicSpecies.') |
|
| 218 |
END IF |
|
| 219 |
END IF |
|
| 220 |
|
|
| 221 |
IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESFORMAT",Search,Clean=".-_")) THEN |
|
| 222 |
Line=Adjustl(Search%Line) |
|
| 223 |
Idx=Index(Line,' ') |
|
| 224 |
Line=AdjustL(Line(Idx+1:)) |
|
| 225 |
Call UpCase(Line) |
|
| 226 |
SELECT CASE (Line) |
|
| 227 |
CASE ('ANG','NOTSCALEDCARTESIANANG')
|
|
| 228 |
Siesta_Unit_Read=1.d0 |
|
| 229 |
CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
|
|
| 230 |
CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
|
|
| 231 |
Siesta_Unit_Read=a0 |
|
| 232 |
END SELECT |
|
| 233 |
END IF |
|
| 234 |
|
|
| 235 |
IF (SearchInput(Siesta_Input,"ATOMICCOORFORMATOUT",Search,Clean=".-_")) THEN |
|
| 236 |
Line=Adjustl(Search%Line) |
|
| 237 |
Idx=Index(Line,' ') |
|
| 238 |
Line=AdjustL(Line(Idx+1:)) |
|
| 239 |
Call UpCase(Line) |
|
| 240 |
SELECT CASE (Line) |
|
| 241 |
CASE ('ANG','NOTSCALEDCARTESIANANG')
|
|
| 242 |
Siesta_Unit_Write=1.d0 |
|
| 243 |
CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
|
|
| 244 |
CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
|
|
| 245 |
Siesta_Unit_Write=Unit |
|
| 246 |
END SELECT |
|
| 247 |
END IF |
|
| 248 |
|
|
| 249 |
IF (SearchInput(Siesta_Input,"LATTICECONSTANT",Search,Clean=".-_")) THEN |
|
| 250 |
Line=Adjustl(Search%Line) |
|
| 251 |
! We discar the label |
|
| 252 |
Idx=Index(Line,' ') |
|
| 253 |
Line=AdjustL(Line(Idx+1:)) |
|
| 254 |
! we read the value |
|
| 255 |
Read(Line,*) Siesta_LatticeConstant |
|
| 256 |
! We discard the value |
|
| 257 |
Idx=Index(Line,' ') |
|
| 258 |
Line=AdjustL(Line(Idx+1:)) |
|
| 259 |
! We read the unit |
|
| 260 |
Call UpCase(Line) |
|
| 261 |
SELECT CASE (Line) |
|
| 262 |
CASE ('ANG')
|
|
| 263 |
Siesta_Lat_Unit=1.d0 |
|
| 264 |
CASE ('BOHR')
|
|
| 265 |
Siesta_Lat_Unit=Unit |
|
| 266 |
END SELECT |
|
| 267 |
! for now! |
|
| 268 |
Call Die('ReadInput_siesta:LatticeConstant',"For now, periodic calculations are NOT possible in Opt'n Path")
|
|
| 269 |
|
|
| 270 |
END IF |
|
| 271 |
|
|
| 272 |
IF (SearchInput(Siesta_Input,"LATTICEVECTORS",Search,Clean=".-_")) THEN |
|
| 273 |
|
|
| 274 |
Call Die('ReadInput_siesta:LatticeVectors',"For now, periodic calculations are NOT possible in Opt'n Path")
|
|
| 275 |
|
|
| 276 |
END IF |
|
| 277 |
|
|
| 278 |
IF (SearchInput(Siesta_Input,"LATTICEPARAMETERS",Search,Clean=".-_")) THEN |
|
| 279 |
|
|
| 280 |
Call Die('ReadInput_siesta:LatticeParameters',"For now, periodic calculations are NOT possible in Opt'n Path")
|
|
| 281 |
|
|
| 282 |
END IF |
|
| 283 |
|
|
| 284 |
IF (SearchInput(Siesta_Input,"SUPERCELL",Search,Clean=".-_")) THEN |
|
| 285 |
|
|
| 286 |
Call Die('ReadInput_siesta:SuperCell',"SuperCell NOT possible in Opt'n Path")
|
|
| 287 |
|
|
| 288 |
END IF |
|
| 289 |
|
|
| 290 |
Call Die('Readinput_siesta:fin',' Lecture finie')
|
|
| 291 |
|
|
| 292 |
if (debug) Call Header("Exiting ReadInput_Siesta")
|
|
| 293 |
|
|
| 294 |
END SUBROUTINE READINPUT_SIESTA |
|
| 69 |
if (debug) Call Header("Exiting WriteInput_Siesta")
|
|
| 70 |
|
|
| 71 |
END SUBROUTINE WRITEINPUT_SIESTA |
|
| src/Extrapol_mixed.f90 (revision 7) | ||
|---|---|---|
| 122 | 122 |
if (s>=dist) THEN |
| 123 | 123 |
s=s-dist |
| 124 | 124 |
IdxGeom=IdxGeom+1 |
| 125 |
SGeom(IdxGeom)=s+IdxGeom*dist |
|
| 125 | 126 |
XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) |
| 126 | 127 |
|
| 127 | 128 |
IntCoordF(IdxGeom,:)=IntCoordTmp |
| ... | ... | |
| 154 | 155 |
|
| 155 | 156 |
END IF |
| 156 | 157 |
END IF |
| 157 |
ENDDO |
|
| 158 |
END DO
|
|
| 158 | 159 |
|
| 159 | 160 |
|
| 160 | 161 |
if ((s>=0.1*dist).AND.(IdxGeom.EQ.NGeomF)) THEN |
| ... | ... | |
| 168 | 169 |
IdxGeom=NGeomF |
| 169 | 170 |
|
| 170 | 171 |
! We have to add the last geometry. We copy the last geom of Initial geometries. |
| 171 |
IntCoordF(IdxGeom,:)=IntCoordI(NGeomI,:)
|
|
| 172 |
IntTangent(IdxGeom,:)=DerInt
|
|
| 172 |
IntCoordF(IdxGeom,:)=IntCoordI(NGeomI,:) |
|
| 173 |
IntTangent(IdxGeom,:)=DerInt |
|
| 173 | 174 |
|
| 174 | 175 |
! we convert it to cartesian geom |
| 175 | 176 |
call Mixed2Cart(Nat,IndZmat,IntCoordF(IdxGeom,1),XyzTmp2(1,1)) |
| src/ReadAnaList.f90 (revision 7) | ||
|---|---|---|
| 1 |
SUBROUTINE ReadAnaList |
|
| 2 |
! This routines read a list of geometrical variables to monitor |
|
| 3 |
! This is inspired from Xyz2Path (that was inspired by Xyz2scan ...) |
|
| 4 |
|
|
| 5 |
|
|
| 6 |
use VarTypes |
|
| 7 |
use Path_module |
|
| 8 |
use Io_module |
|
| 9 |
|
|
| 10 |
IMPLICIT NONE |
|
| 11 |
|
|
| 12 |
|
|
| 13 |
INTERFACE |
|
| 14 |
function valid(string) result (isValid) |
|
| 15 |
CHARACTER(*), intent(in) :: string |
|
| 16 |
logical :: isValid |
|
| 17 |
END function VALID |
|
| 18 |
|
|
| 19 |
SUBROUTINE die(routine, msg, file, line, unit) |
|
| 20 |
|
|
| 21 |
Use VarTypes |
|
| 22 |
Use io_module |
|
| 23 |
|
|
| 24 |
implicit none |
|
| 25 |
|
|
| 26 |
character(len=*), intent(in) :: routine, msg |
|
| 27 |
character(len=*), intent(in), optional :: file |
|
| 28 |
integer(KINT), intent(in), optional :: line, unit |
|
| 29 |
|
|
| 30 |
END SUBROUTINE die |
|
| 31 |
|
|
| 32 |
END INTERFACE |
|
| 33 |
|
|
| 34 |
LOGICAL :: Debug,FirstVar,FirstCom |
|
| 35 |
INTEGER(KINT) :: At1,At2,At3,At4 |
|
| 36 |
INTEGER(KINT) :: I,Idx,J |
|
| 37 |
|
|
| 38 |
CHARACTER(LCHARS) :: Line |
|
| 39 |
|
|
| 40 |
Debug=Valid('ReadAnaList')
|
|
| 41 |
If (Debug) Call Header("Entering ReadAnaList")
|
|
| 42 |
|
|
| 43 |
Allocate(Bary) |
|
| 44 |
Nullify(Bary%Next) |
|
| 45 |
CurBary => Bary |
|
| 46 |
|
|
| 47 |
ALLOCATE(GeomList) |
|
| 48 |
Nullify(GeomList%Next) |
|
| 49 |
CurVar => GeomList |
|
| 50 |
|
|
| 51 |
NbCom=0 |
|
| 52 |
NbVar=0 |
|
| 53 |
FormAna='(1X' |
|
| 54 |
DO I=1, Nb |
|
| 55 |
READ(IOIN,'(A)') Line |
|
| 56 |
Line=AdjustL(Line) |
|
| 57 |
Call Upcase(Line) |
|
| 58 |
|
|
| 59 |
SELECT CASE (Line(1:1)) |
|
| 60 |
CASE ('B')
|
|
| 61 |
! this is a bond |
|
| 62 |
CurVar%Type="BOND" |
|
| 63 |
Idx=Index(Line," ") |
|
| 64 |
Line=Line(Idx+1:) |
|
| 65 |
Read(Line,*) At1, At2 |
|
| 66 |
CurVar%At1=At1 |
|
| 67 |
CurVar%At2=At2 |
|
| 68 |
CurVar%Value=0. |
|
| 69 |
! CurVar%PrintFactor=1. |
|
| 70 |
Curvar%SignDihedral=1 |
|
| 71 |
if (debug) THEN |
|
| 72 |
WRITE(*,'("# b ",2I3)') At1,At2
|
|
| 73 |
END IF |
|
| 74 |
FormAna=TRIM(FormAna) //',1X,F7.3' |
|
| 75 |
Allocate(CurVar%Next) |
|
| 76 |
CurVar => CurVar%Next |
|
| 77 |
Nullify(CurVar%Next) |
|
| 78 |
CASE ('A')
|
|
| 79 |
! this is a valence angle |
|
| 80 |
CurVar%Type="ANGLE" |
|
| 81 |
Idx=Index(Line," ") |
|
| 82 |
Line=Line(Idx+1:) |
|
| 83 |
Read(Line,*) At1, At2,At3 |
|
| 84 |
CurVar%At1=At1 |
|
| 85 |
CurVar%At2=At2 |
|
| 86 |
CurVar%At3=At3 |
|
| 87 |
CurVar%Value=0. |
|
| 88 |
! CurVar%PrintFactor=180./Pi |
|
| 89 |
Curvar%SignDihedral=1 |
|
| 90 |
if (debug) THEN |
|
| 91 |
WRITE(*,'("# a ",4(I3))') At1,At2,At3
|
|
| 92 |
END IF |
|
| 93 |
FormAna=TRIM(FormAna) //',1X,F7.2' |
|
| 94 |
Allocate(CurVar%Next) |
|
| 95 |
CurVar => CurVar%Next |
|
| 96 |
Nullify(CurVar%Next) |
|
| 97 |
CASE ('D')
|
|
| 98 |
! this is a dihedral |
|
| 99 |
CurVar%Type="DIHEDRAL" |
|
| 100 |
Idx=Index(Line," ") |
|
| 101 |
Line=Line(Idx+1:) |
|
| 102 |
Read(Line,*) At1, At2,At3,At4 |
|
| 103 |
CurVar%At1=At1 |
|
| 104 |
CurVar%At2=At2 |
|
| 105 |
CurVar%At3=At3 |
|
| 106 |
CurVar%At4=At4 |
|
| 107 |
CurVar%Value=0. |
|
| 108 |
! CurVar%PrintFactor=180./Pi |
|
| 109 |
Curvar%SignDihedral=1 |
|
| 110 |
if (debug) THEN |
|
| 111 |
WRITE(*,'("# d ",4(I3))') At1,At2,At3,At4
|
|
| 112 |
END IF |
|
| 113 |
Allocate(CurVar%Next) |
|
| 114 |
CurVar => CurVar%Next |
|
| 115 |
Nullify(CurVar%Next) |
|
| 116 |
CASE ('C')
|
|
| 117 |
NbCom=NbCom+1 |
|
| 118 |
Idx=Index(Line," ") |
|
| 119 |
Line=Line(Idx+1:) |
|
| 120 |
READ(Line,*) At1 |
|
| 121 |
Allocate(CurBary%ListAtoms(0:At1)) |
|
| 122 |
Allocate(CurBary%Weights(1:At1)) |
|
| 123 |
CurBary%Weights=1. |
|
| 124 |
CurBary%ListAtoms(0)=At1 |
|
| 125 |
Idx=Index(Line," ") |
|
| 126 |
Line=Line(Idx+1:) |
|
| 127 |
READ(Line,*) CurBary%ListAtoms(1:At1) |
|
| 128 |
if (debug) THEN |
|
| 129 |
WRITE(*,'("# c ",I4,20(I3))') At1, &
|
|
| 130 |
(CurBary%ListAtoms(j),j=1,At1) |
|
| 131 |
END IF |
|
| 132 |
FormAna=Trim(FormAna) //',1X,F7.2' |
|
| 133 |
|
|
| 134 |
Allocate(CurBary%Next) |
|
| 135 |
CurBary => CurBary%Next |
|
| 136 |
Nullify(CurBary%Next) |
|
| 137 |
|
|
| 138 |
CASE Default |
|
| 139 |
Call Die('ReadAnaList','Variable not recognized: ' //Line, Unit=IOOUT)
|
|
| 140 |
END SELECT |
|
| 141 |
|
|
| 142 |
END DO |
|
| 143 |
|
|
| 144 |
NbVar=Nb-NbCom |
|
| 145 |
|
|
| 146 |
CurVar => GeomList |
|
| 147 |
I=0 |
|
| 148 |
|
|
| 149 |
Allocate(PrintGeomFactor(NbVar)) |
|
| 150 |
DO WHILE (associated(CurVar%Next)) |
|
| 151 |
I=I+1 |
|
| 152 |
SELECT CASE (CurVar%TYPE) |
|
| 153 |
CASE ('ANGLE','DIHEDRAL')
|
|
| 154 |
PrintGeomFactor(I)=180./pi |
|
| 155 |
CASE DEFAULT |
|
| 156 |
PrintGeomFactor(I)=1. |
|
| 157 |
END SELECT |
|
| 158 |
CurVar => CurVar%Next |
|
| 159 |
END DO |
|
| 160 |
|
|
| 161 |
FormAna=Trim(FormAna) // ',1X,F7.2,1X,F15.6)' |
|
| 162 |
|
|
| 163 |
|
|
| 164 |
! if (debug) WRITE(*,*) "NbVar, NbCom,Nb",NbVar,NbCom,Nb |
|
| 165 |
|
|
| 166 |
if (debug) Call Header(" Exiting ReadAnaList")
|
|
| 167 |
|
|
| 168 |
END SUBROUTINE ReadAnaList |
|
| src/Path.f90 (revision 7) | ||
|---|---|---|
| 212 | 212 |
! We add Siesta (and start thinking of add CP2K) |
| 213 | 213 |
! as a new "energy engine" code |
| 214 | 214 |
! TO DO to go to 1.5: finish cleaning... |
| 215 |
! |
|
| 215 | 216 |
! 2013 Feb Opt'n Path v1.48 |
| 216 | 217 |
! We add some keyword for more output for Geometry Optimizations |
| 217 | 218 |
! GeomFile: name of the file to print out the geometries and their energy |
| 218 | 219 |
! as Xyz file. (only format for now) |
| 220 |
! |
|
| 221 |
! 2013 Feb opt'n Path v1.49 |
|
| 222 |
! We add the possibility to analyze the geometries in terms |
|
| 223 |
! of Bonds, Angles, dihedrals |
|
| 224 |
! with the option to add CenterOfMass atoms for the analysis (a la Xyz2Path) |
|
| 225 |
! This is basicaly the merging of AnaPath into Opt'n Path... at last ! |
|
| 226 |
! This is done by adding a Flag into the Path namelist: |
|
| 227 |
! AnaGeom: logical, if true Path looks for the AnaList namelist |
|
| 228 |
! AnaList: list of variables, if present and if AnaGeom=T then Opt'n Path |
|
| 229 |
! will read it and print the values of the variable in a Path.dat file |
|
| 230 |
! If AnaGeom is T but Analist is not given, then Path.dat just contains the energy. |
|
| 219 | 231 |
|
| 220 | 232 |
PROGRAM PathOpt |
| 221 | 233 |
|
| ... | ... | |
| 249 | 261 |
LOGICAL :: FChkFrozen |
| 250 | 262 |
|
| 251 | 263 |
! Energies for all points along the path at the previous iteration |
| 252 |
REAL(KREAL), ALLOCATABLE :: EPathold(:) ! NGeomF, where it is deallocated: Prakash
|
|
| 264 |
REAL(KREAL), ALLOCATABLE :: EPathold(:) ! NGeomF |
|
| 253 | 265 |
! Maximum step allowed for this geometry |
| 254 |
REAL(KREAL), ALLOCATABLE :: MaxStep(:) ! NGeomF, where it is deallocated: Prakash
|
|
| 266 |
REAL(KREAL), ALLOCATABLE :: MaxStep(:) ! NGeomF |
|
| 255 | 267 |
|
| 256 | 268 |
! these are used to read temporary coordinates |
| 257 | 269 |
LOGICAL :: FFrozen,FCart |
| ... | ... | |
| 272 | 284 |
|
| 273 | 285 |
! INTEGER(KINT), EXTERNAL :: Iargc |
| 274 | 286 |
INTEGER(KINT), EXTERNAL :: ConvertNumAt |
| 287 |
INTEGER(KINT) :: IoS |
|
| 275 | 288 |
|
| 276 | 289 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| 277 | 290 |
! |
| ... | ... | |
| 327 | 340 |
OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, & |
| 328 | 341 |
FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, & |
| 329 | 342 |
NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, & |
| 330 |
NGintMax, Align, CalcEReac,CalcEProd, GeomFile |
|
| 343 |
NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile
|
|
| 331 | 344 |
|
| 332 | 345 |
Namelist/cartlist/list |
| 333 | 346 |
Namelist/frozenlist/list |
| 347 |
Namelist/analist/nb |
|
| 334 | 348 |
|
| 335 | 349 |
|
| 336 | 350 |
Flag_Opt_Geom = .FALSE. ! Initialized. |
| ... | ... | |
| 368 | 382 |
WRITE(*,*) "Input: string that indicates the type of the input geometries." |
| 369 | 383 |
WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp" |
| 370 | 384 |
WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations." |
| 371 |
WRITE(*,*) " Accepted values are: Gaussian, Vasp, Mopac, Test, Siesta or Ext"
|
|
| 372 |
WRITE(*,*) " In case of a Gaussian or Mopac calculations, input must be set to Cart." |
|
| 385 |
WRITE(*,*) " Accepted values are: Gaussian, Vasp, Mopac, TurboMole, Siesta, Test or Ext"
|
|
| 386 |
WRITE(*,*) " * In case of a Gaussian or Mopac calculations, input must be set to Cart."
|
|
| 373 | 387 |
WRITE(*,*) " One example of a gaussian/mopac input should be added at the end of the" |
| 374 |
WRITE(*,*) " input file.See example file Test_G03.path or Test_Mopac.path"
|
|
| 375 |
WRITE(*,*) " In the case of a VASP calculation, if input is set to Cart, then" |
|
| 376 |
WRITE(*,*) " the preamble of a VASP calculation should be added at the end of" |
|
| 377 |
WRITE(*,*) " the input file. See example file Test_VASP_cart.path" |
|
| 378 |
WRITE(*,*) " In the case of a VASP calculation, one should also give value of the " |
|
| 388 |
WRITE(*,*) " input file. See example file Test_G03.path or Test_Mopac.path"
|
|
| 389 |
WRITE(*,*) " * In the case of a VASP calculation, if input is set to Cart, then"
|
|
| 390 |
WRITE(*,*) " the preamble of a VASP calculation should be added at the end of"
|
|
| 391 |
WRITE(*,*) " the input file. See example file Test_VASP_cart.path"
|
|
| 392 |
WRITE(*,*) " In the case of a VASP calculation, one should also give value of the "
|
|
| 379 | 393 |
WRITE(*,*) " Runmode variable" |
| 394 |
WRITE(*,*) " * In the case of a SIESTA calculation, an example of a Siesta input file" |
|
| 395 |
WRITE(*,*) " should be added at the end of the input file. See Test_Siesta.path" |
|
| 380 | 396 |
WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy" |
| 381 | 397 |
WRITE(*,*) " and gradient of the whole path or not. If one wants to use VASP," |
| 382 | 398 |
WRITE(*,*) " Runmode must be set to PARA." |
| ... | ... | |
| 413 | 429 |
WRITE(*,*) "Smax: Maximum length of a step during path optimization" |
| 414 | 430 |
WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary" |
| 415 | 431 |
WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken" |
| 416 |
WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path"
|
|
| 432 |
WRITE(*,*) "FTan: When moving the path, this gives the proportion of the displacement tangent to the path"
|
|
| 417 | 433 |
WRITE(*,*) " that is kept. FTan=1. corresponds to the full displacement, " |
| 418 |
WRITE(*,*) " whereas FTan=0. gives a displacement orthogonal to the path." |
|
| 434 |
WRITE(*,*) " whereas FTan=0. gives a displacement orthogonal to the path. Default: Ftan=0."
|
|
| 419 | 435 |
WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization." |
| 420 | 436 |
WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path." |
| 421 | 437 |
WRITE(*,*) " This option indicates the first step where spline interpolation is used." |
| ... | ... | |
| 451 | 467 |
WRITE(*,*) " If True, a &cartlist namelist containing the list of cart atoms must be given." |
| 452 | 468 |
WRITE(*,*) " By default, only frozen atoms are described in cartesian coordinates." |
| 453 | 469 |
WRITE(*,*) "" |
| 470 |
WRITE(*,*) "AnaGEom: If TRUE, Opt'n Path will create a file .dat for geometries analysis." |
|
| 471 |
WRITE(*,*) " If True, Opt'n Path will look for the AnaList namelist after the Path Namelist." |
|
| 472 |
WRITE(*,*) "AnaList: list of variables for geometry analysis. If present and if AnaGeom=T then" |
|
| 473 |
WRITE(*,*) " Opt'n Path will read it and print the values of the variable in a .dat file" |
|
| 474 |
WRITE(*,*) " If AnaGeom is T but Analist is not given, then Path.dat just contains the energy." |
|
| 475 |
WRITE(*,*) "Analist contains the number of variables to monitor: Nb, and is followed by the description" |
|
| 476 |
WRITE(*,*) "of the variables among:" |
|
| 477 |
WRITE(*,*) "b(ond) At1 At2" |
|
| 478 |
WRITE(*,*) "a(ngle) At1 At2 At3" |
|
| 479 |
WRITE(*,*) "d(ihedral) At1 At2 At3 At4" |
|
| 480 |
WRITE(*,*) "c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added" |
|
| 481 |
WRITE(*,*) "at the end of the real atoms of the system" |
|
| 482 |
WRITE(*,*) "" |
|
| 454 | 483 |
WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry." |
| 455 | 484 |
WRITE(*,*) " If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. " |
| 456 | 485 |
WRITE(*,*) " It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]" |
| ... | ... | |
| 470 | 499 |
END SELECT |
| 471 | 500 |
|
| 472 | 501 |
|
| 473 |
|
|
| 474 | 502 |
! We initiliaze variables |
| 475 | 503 |
Pi=dacos(-1.0d0) |
| 476 | 504 |
Nat=1 |
| ... | ... | |
| 507 | 535 |
CalcEReac=.FALSE. |
| 508 | 536 |
CalcEProd=.FALSE. |
| 509 | 537 |
EReac=0.d0 |
| 510 |
EProd=0.d0 |
|
| 538 |
EProd=0.d0
|
|
| 511 | 539 |
OptGeom=-1 |
| 512 | 540 |
GeomFile="EMPTY" |
| 541 |
AnaGeom=.FALSE. |
|
| 542 |
AnaFile="EMPTY" |
|
| 543 |
Nb=0 |
|
| 544 |
NbCom=0 |
|
| 513 | 545 |
PathOnly=.FALSE. |
| 514 | 546 |
Prog='EMPTY' |
| 515 | 547 |
ProgExe='EMPTY' |
| ... | ... | |
| 580 | 612 |
END IF |
| 581 | 613 |
END IF |
| 582 | 614 |
|
| 583 |
If (COORD.EQ.'XYZ') THEN |
|
| 615 |
If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN |
|
Formats disponibles : Unified diff