root / src / egrad_test_2D.f90 @ 5
History  View  Annotate  Download (4.4 kB)
1 
SUBROUTINE egrad_test_2D(nat,e,geom,grad) 

2 
! This program computes the energy and the gradient 
3 
! to cheat PATH we use a 3 atoms system 
4  
5 
IMPLICIT NONE 
6 
integer, parameter :: KINT = kind(1) 
7 
integer, parameter :: KREAL = kind(1.0d0) 
8  
9 
INTEGER(KINT), INTENT(IN) :: Nat 
10 
REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3) 
11 
REAL(KREAL), INTENT(IN) :: Geom(Nat,3) 
12  
13 
! Bohr > Angstr 
14 
! real(KREAL), parameter :: BOHR = 0.52917726D+00 
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.01d0 
21 
INTEGER(KINT), PARAMETER :: IOTMP=99 
22 
! Variables 
23 
INTEGER(KINT) :: i, iat, jat,j 
24 
REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:) 
25  
26 
LOGICAL :: Debug 
27 
LOGICAL, save :: First 
28 

29 
REAL(KREAL) :: x,y 
30  
31  
32  
33 
INTERFACE 
34 
function valid(string) result (isValid) 
35 
CHARACTER(*), intent(in) :: string 
36 
logical :: isValid 
37 
END function VALID 
38  
39  
40 
function E2D(natoms,Xyz) 
41  
42 
use Path_module, only : order 
43 
use Io_module, only : IOTMP 
44  
45 
IMPLICIT NONE 
46 
integer, parameter :: KINT = kind(1) 
47 
integer, parameter :: KREAL = kind(1.0d0) 
48  
49  
50 
INTEGER(KINT) ,INTENT(IN) :: natoms 
51 
REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3) 
52 
REAL(KREAL) :: E2D 
53 
END function E2D 
54  
55 

56 
END INTERFACE 
57  
58 
debug=valid('egrad_test_2D') 
59  
60 
if (debug) WRITE(*,*) '================ Entering Egrad_test ===================' 
61 
if (debug) THEN 
62 
WRITE(*,*) "Cartesian Geometry" 
63 
DO I=1,Nat 
64 
WRITE(*,'(1X,I5,3(1X,F12.6))') I,Geom(i,1:3) 
65 
END DO 
66 
END IF 
67 

68 
ALLOCATE(XyZTmp(Nat,3),GradTmp(Nat,3)) 
69  
70 
if (first) THEN 
71 
OPEN(IOTMP,File="tt") 
72 
DO I=1, 501 
73 
DO j=1,501 
74 
x=(i1.)/100. 
75 
y=(j1.)/100. 
76 
xyztmp(2,1)=x 
77 
xyztmp(3,1)=y 
78 
write(IOTMP,*) x,y,E2D(nat,xyztmp) 
79 
END DO 
80 
! we skip a blank line 
81 
WRITe(IOTMP,*) 
82 
END DO 
83 
CLOSE(IOTMP) 
84 
First=.FALSE. 
85 
END IF 
86  
87 
e=E2D(nat,Geom) 
88  
89 
! We now calculate the gradients using numerical derivatives 
90 
Grad=0.d0 
91 
GradTmp=0.d0 
92 
do iat=1,1 
93 
do jat=2,nat 
94 
xyztmp=geom 
95 
do i=1,6 
96 
xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh 
97 
gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*E2D(nat,xyztmp) 
98 
end do 
99 
end do 
100 
end do 
101 
gradTmp=gradTmp/(60.*hh) 
102 
do iat=1,nat 
103 
grad(3*iat2:3*iat)=gradTmp(iat,1:3) 
104 
end do 
105  
106 
if (debug) THEN 
107 
WRITE(*,*) "Cartesian gradient GradTmp" 
108 
DO I=1,Nat 
109 
WRITE(*,'(1X,I5,3(1X,F12.6))') I,Gradtmp(i,1:3) 
110 
END DO 
111 
WRITE(*,*) "Cartesian gradient Grad" 
112 
WRITE(*,'(50(1X,F12.6))') Grad 
113  
114 
END IF 
115  
116  
117 
deallocate(xyztmp,gradTmp) 
118  
119  
120 
if (debug) WRITE(*,*) '================ Exiting Egrad_test_2D ===================' 
121  
122 
! ====================================================================== 
123 
end SUBROUTINE egrad_test_2D 
124  
125  
126  
127 
function E2D(natoms,Xyz) 
128  
129 
use Path_module, only : order 
130 
use Io_module, only : IOTMP 
131  
132 
IMPLICIT NONE 
133 
integer, parameter :: KINT = kind(1) 
134 
integer, parameter :: KREAL = kind(1.0d0) 
135  
136  
137 
INTEGER(KINT) ,INTENT(IN) :: natoms 
138 
REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3) 
139 
REAL(KREAL) :: E,x,y, E2D 
140 
REAL(KREAL), save :: Data(501,501) 
141  
142 
INTEGER(KINT) :: i,j 
143 
REAL(KREAL) :: dCH,dNH,dCN 
144  
145 
! real(KREAL) :: autoA,autoeV 
146 
real(KREAL) :: autoeV 
147 
! parameter (autoA=0.52917715d0) 
148 
parameter (autoeV=27.21183d0) 
149  
150  
151 
real(KREAL) :: r1, r2 
152  
153 
LOGICAL :: Debug 
154 
LOGICAL, SAVE :: First=.TRUE. 
155 

156 
INTERFACE 
157 
function valid(string) result (isValid) 
158 
CHARACTER(*), intent(in) :: string 
159 
logical :: isValid 
160 
END function VALID 
161 

162 
END INTERFACE 
163  
164 
debug=valid('e2D') 
165  
166 
if (First) THEN 
167 
! we read the data 
168 
open(IOTMP,File="CCSDT.subtr.reg") 
169 
DO I=1, 501 
170 
DO j=1,501 
171 
READ(IOTMP,*) x,y,E 
172 
Data(int(x*100+0.1)+1,int(y*100+0.1)+1)=E 
173 
write(*,*) x,y,int(x*100+0.1)+1, int(y*100+0.1)+1,E 
174 
END DO 
175 
! we skip a blank line 
176 
READ(IOTMP,*) 
177 
END DO 
178 
CLOSE(IOTMP) 
179 
First=.FALSE. 
180 
END IF 
181  
182 
x=Xyz(2,1) 
183 
y=Xyz(3,1) 
184 
E2D=Data(int(x*100+0.1)+1,int(y*100+0.1)+1) 
185  
186 
return 
187 
end function E2D 