root / src / egrad_LEPS.f90 @ 8
History  View  Annotate  Download (3.8 kB)
1 
SUBROUTINE egrad_LEPS(nat,e,geom,grad) 

2 
! This program computes the ernergy and gradien in cartesian coordinates 
3 
! for the cartesian geometry Geom 
4  
5 
IMPLICIT NONE 
6 
integer, parameter :: KINT = kind(1) 
7 
integer, parameter :: KREAL = kind(1.0d0) 
8  
9 
! Number of atoms 
10 
INTEGER(KINT), INTENT(IN) :: Nat 
11 
! Input geometry 
12 
REAL(KREAL), INTENT(IN) :: Geom(Nat,3) 
13 
! output energy and gradient 
14 
REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3) 
15  
16  
17 
! Parameters to define the surface 
18 
INTEGER(KINT), DIMENSION(6), PARAMETER :: IECOEF = (/1,9,45,45,9,1/) 
19 
INTEGER(KINT), DIMENSION(6), PARAMETER :: ISCOEF = (/3,2,1,1,2,3/) 
20 
REAL(KREAL), PARAMETER :: hh=0.001d0 
21 

22 
! Variables 
23 
INTEGER(KINT) :: i,iat,jat 
24 
REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:) 
25 

26 
LOGICAL :: Debug 
27 

28 
REAL(KREAL), external :: ELEPS_xyz 
29  
30  
31 
INTERFACE 
32 
function valid(string) result (isValid) 
33 
CHARACTER(*), intent(in) :: string 
34 
logical :: isValid 
35 
END function VALID 
36 

37 
END INTERFACE 
38  
39 
debug=valid('egrad_leps') 
40  
41 
if (debug) WRITE(*,*) '================ Entering Egrad_leps ===================' 
42 
if (debug) THEN 
43 
WRITE(*,*) "Cartesian Geometry" 
44 
DO I=1,Nat 
45 
WRITE(*,'(1X,I5,3(1X,F12.6))') I,Geom(i,1:3) 
46 
END DO 
47 
END IF 
48 

49 
ALLOCATE(XyZTmp(Nat,3),GradTmp(Nat,3)) 
50  
51 
e=ELEPS_xyz(nat,Geom) 
52  
53 
! We now calculate the gradients using numerical derivatives 
54 
Grad=0.d0 
55 
GradTmp=0.d0 
56 
do iat=1,3 
57 
do jat=1,nat 
58 
xyztmp=geom 
59 
do i=1,6 
60 
xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh 
61 
gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*ELEPS_xyz(nat,xyztmp) 
62 
end do 
63 
end do 
64 
end do 
65 
gradTmp=gradTmp/(60.*hh) 
66 
do iat=1,nat 
67 
grad(3*iat2:3*iat)=gradTmp(iat,1:3) 
68 
end do 
69  
70 
if (debug) THEN 
71 
WRITE(*,*) "Cartesian gradient " 
72 
DO I=1,Nat 
73 
WRITE(*,'(1X,I5,3(1X,F12.6))') I,Gradtmp(i,1:3) 
74 
END DO 
75 
END IF 
76  
77  
78 
deallocate(xyztmp,gradTmp) 
79  
80  
81 
if (debug) WRITE(*,*) '================ Exiting Egrad_leps ===================' 
82  
83 
! ====================================================================== 
84 
end 
85  
86  
87  
88 
function ELEPS_xyz(natoms,Xyz) 
89  
90 
use Path_module, only : order 
91 
use Io_module, only : au2ev 
92  
93 
IMPLICIT NONE 
94 
integer, parameter :: KINT = kind(1) 
95 
integer, parameter :: KREAL = kind(1.0d0) 
96  
97  
98 
INTEGER(KINT) ,INTENT(IN) :: natoms 
99 
REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3) 
100 
REAL(KREAL) :: ELEPS_xyz 
101  
102 
INTEGER(KINT) :: i 
103 
REAL(KREAL) :: rAB, rBC, rAC 
104  
105 
REAL(KREAL), PARAMETER :: a=0.05, b=0.30, c=0.05 
106 
REAL(KREAL), PARAMETER :: dAB=4.746, dBC=4.747, dAC=3.445 
107 
REAL(KREAL), PARAMETER :: r0=0.742, alpha=1.942 
108  
109 
REAL(KREAL) :: Qbc,QAc,Qab,Jab,Jbc,Jac 
110  
111 
rAB=0. 
112 
rAC=0. 
113 
rBC=0. 
114 
do i=1,3 
115 
rBC=rBC+(xyz(Order(3),i)xyz(order(2),i))**2 
116 
rAB=rAB+(xyz(order(1),i)xyz(order(2),i))**2 
117 
rAC=rAC+(xyz(order(1),i)xyz(order(3),i))**2 
118 
end do 
119 
rAB=sqrt(rAB) 
120 
rBC=sqrt(rBC) 
121 
rAC=sqrt(rAC) 
122  
123 
! For consistency with previous articles : rAC is constrained to be 3.742 
124 
rAC=3.742 
125  
126 
Qab = dab*(1.5*exp(2*alpha*(rABr0))exp(alpha*(rABr0)))/(2+2*a) 
127 
Qbc = dbc*(1.5*exp(2*alpha*(rBCr0))exp(alpha*(rBCr0)))/(2+2*b) 
128 
Qac = dac*(1.5*exp (2*alpha*(rACr0))exp(alpha*(rACr0)))/(2+2*c) 
129 
Jab = dab*(exp(2*alpha*(rABr0))6*exp(alpha*(rABr0)))/(4+4*a) 
130 
Jbc = dbc*(exp(2*alpha*(rBCr0))6*exp(alpha*(rBCr0)))/(4+4*b) 
131 
Jac = dac*(exp(2*alpha*(rACr0))6*exp(alpha*(rACr0)))/(4+4*c) 
132  
133 
! V(x,y) = Qab(x)+Qbc(y)+Qac(rac)sqrt(Jab(x)**2+Jbc(y)**2+Jac(x)**2Jab(x)*Jbc(y)Jbc(y)*Jac(rac)Jab(x)*Jac(rac)) 
134  
135  
136 
ELEPS_xyz=(Qab+Qbc+Qacsqrt(Jab**2+Jbc**2+Jac**2Jab*JbcJbc*JacJab*Jac))/au2eV 
137  
138 
return 
139 
end function ELEPS_xyz 