Line data Source code
1 : !!****m* ABINIT/m_pred_srkhna14
2 : !! NAME
3 : !! m_pred_srkna14
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, JCC, SE)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_pred_srkhna14
22 :
23 : use defs_basis
24 : use m_abicore
25 : use m_abimover
26 : use m_abihist
27 :
28 : use m_geometry, only : xcart2xred, xred2xcart, metric
29 :
30 : implicit none
31 :
32 : private
33 : !!***
34 :
35 : public :: pred_srkna14
36 : !!***
37 :
38 : contains
39 : !!***
40 :
41 : !!****f* ABINIT/pred_srkna14
42 : !! NAME
43 : !! pred_srkna14
44 : !!
45 : !! FUNCTION
46 : !! Ionmov predictors (14) Srkna14 molecular dynamics
47 : !!
48 : !! IONMOV 14:
49 : !! Simple molecular dynamics with a symplectic algorithm proposed
50 : !! by S.Blanes and P.C.Moans, called SRKNa14 in Practical symplectic partitioned
51 : !! Runge--Kutta and Runge--Kutta--Nystrom methods, Journal of Computational
52 : !! and Applied Mathematics archive, volume 142, issue 2 (May 2002), pages 313 - 330 [[cite:Blanes2002]].
53 : !! of the kind first published by H. Yoshida, Construction of higher order symplectic
54 : !! integrators, Physics Letters A, volume 150, number 5 to 7, pages 262 - 268 [[cite:Yoshida1990]]
55 : !! This algorithm requires at least 14 evaluation of the forces (actually 15 are done
56 : !! within Abinit) per time step. At this cost it usually gives much better
57 : !! energy conservation than the verlet algorithm (ionmov 6) for a 30 times bigger
58 : !! value of <a href="varrlx.html#dtion">dtion</a>. Notice that the potential
59 : !! energy of the initial atomic configuration is never evaluated using this
60 : !! algorithm.
61 : !!
62 : !! INPUTS
63 : !! ab_mover <type(abimover)> : Datatype with all the information needed by the preditor
64 : !! itime : Index of the present iteration
65 : !! ntime : Maximal number of iterations
66 : !! icycle : Index of the present cycle
67 : !! zDEBUG : if true print some debugging information
68 : !!
69 : !! SIDE EFFECTS
70 : !! hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
71 : !!
72 : !! SOURCE
73 :
74 16 : subroutine pred_srkna14(ab_mover,hist,icycle,zDEBUG,iexit,skipcycle)
75 :
76 : !Arguments ------------------------------------
77 : !scalars
78 : type(abimover),intent(in) :: ab_mover
79 : type(abihist),intent(inout) :: hist
80 : integer,intent(in) :: icycle
81 : integer,intent(in) :: iexit
82 : logical,intent(in) :: zDEBUG
83 : logical,intent(out) :: skipcycle
84 :
85 : !Local variables-------------------------------
86 : !scalars
87 : integer :: ihist_prev,ii,jj,kk
88 : real(dp) :: ucvol,ucvol_next
89 : real(dp),parameter :: v2tol=tol8
90 : real(dp) :: etotal
91 : logical :: jump_end_of_cycle=.FALSE.
92 : ! character(len=5000) :: message
93 : !arrays
94 : real(dp),save :: aa(15),bb(15)
95 : real(dp) :: acell(3),acell_next(3)
96 : real(dp) :: rprimd(3,3),rprimd_next(3,3)
97 : real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3)
98 32 : real(dp) :: fcart(3,ab_mover%natom),fcart_m(3,ab_mover%natom)
99 32 : real(dp) :: xcart(3,ab_mover%natom)
100 32 : real(dp) :: xred(3,ab_mover%natom)
101 1 : real(dp) :: vel(3,ab_mover%natom)
102 : real(dp) :: strten(6)
103 :
104 : !***************************************************************************
105 : !Beginning of executable session
106 : !***************************************************************************
107 :
108 16 : if(iexit/=0)then
109 : return
110 : end if
111 :
112 15 : jump_end_of_cycle=.FALSE.
113 255 : fcart_m(:,:)=zero
114 :
115 : !write(std_out,*) 'srkna14 03',jump_end_of_cycle
116 : !##########################################################
117 : !### 03. Obtain the present values from the history
118 :
119 15 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
120 15 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
121 :
122 15 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
123 :
124 255 : fcart(:,:)=hist%fcart(:,:,hist%ihist)
125 105 : strten(:) =hist%strten(:,hist%ihist)
126 255 : vel(:,:) =hist%vel(:,:,hist%ihist)
127 15 : etotal =hist%etot(hist%ihist)
128 :
129 15 : if(zDEBUG)then
130 0 : write (std_out,*) 'fcart:'
131 0 : do kk=1,ab_mover%natom
132 0 : write (std_out,*) fcart(:,kk)
133 : end do
134 0 : write (std_out,*) 'vel:'
135 0 : do kk=1,ab_mover%natom
136 0 : write (std_out,*) vel(:,kk)
137 : end do
138 0 : write (std_out,*) 'strten:'
139 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
140 0 : write (std_out,*) 'etotal:'
141 0 : write (std_out,*) etotal
142 : end if
143 :
144 15 : write(std_out,*) 'RMET'
145 60 : do ii=1,3
146 60 : write(std_out,*) rmet(ii,:)
147 : end do
148 :
149 : !write(std_out,*) 'srkna14 04',jump_end_of_cycle
150 : !##########################################################
151 : !### 04. Compute the next values (Only for the first cycle)
152 :
153 15 : if (icycle==1) then
154 :
155 1 : if(zDEBUG) then
156 0 : write(std_out,*) 'Entering only for first cycle'
157 : end if
158 :
159 1 : aa(1) = 0.0378593198406116_dp;
160 1 : aa(2) = 0.102635633102435_dp;
161 1 : aa(3) = -0.0258678882665587_dp;
162 1 : aa(4) = 0.314241403071447_dp;
163 1 : aa(5) = -0.130144459517415_dp;
164 1 : aa(6) = 0.106417700369543_dp;
165 1 : aa(7) = -0.00879424312851058_dp;
166 : aa(8) = 1._dp -&
167 1 : & 2._dp*(aa(1)+aa(2)+aa(3)+aa(4)+aa(5)+aa(6)+aa(7));
168 1 : aa(9) = aa(7);
169 1 : aa(10)= aa(6);
170 1 : aa(11)= aa(5);
171 1 : aa(12)= aa(4);
172 1 : aa(13)= aa(3);
173 1 : aa(14)= aa(2);
174 1 : aa(15)= aa(1);
175 :
176 1 : bb(1) = 0.0_dp
177 1 : bb(2) = 0.09171915262446165_dp;
178 1 : bb(3) = 0.183983170005006_dp;
179 1 : bb(4) = -0.05653436583288827_dp;
180 1 : bb(5) = 0.004914688774712854_dp;
181 1 : bb(6) = 0.143761127168358_dp;
182 1 : bb(7) = 0.328567693746804_dp;
183 1 : bb(8) = 0.5_dp - (bb(1)+bb(2)+bb(3)+bb(4)+bb(5)+bb(6)+bb(7));
184 1 : bb(9) = 0.5_dp - (bb(1)+bb(2)+bb(3)+bb(4)+bb(5)+bb(6)+bb(7));
185 1 : bb(10)= bb(7);
186 1 : bb(11)= bb(6);
187 1 : bb(12)= bb(5);
188 1 : bb(13)= bb(4);
189 1 : bb(14)= bb(3);
190 1 : bb(15)= bb(2);
191 :
192 : acell_next(:)=acell(:)
193 1 : ucvol_next=ucvol
194 : rprimd_next(:,:)=rprimd(:,:)
195 :
196 : ! step 1 of 15
197 :
198 : ! Convert input xred (reduced coordinates) to xcart (cartesian)
199 1 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
200 :
201 17 : vel(:,:) = vel(:,:) + bb(1) * ab_mover%dtion * fcart_m(:,:)
202 :
203 4 : do ii=1,3
204 16 : do jj=1,ab_mover%natom
205 12 : write(std_out,*) xcart(ii,jj), ab_mover%dtion, aa(1), vel(ii,jj)
206 12 : xcart(ii,jj) = xcart(ii,jj) + ab_mover%dtion * aa(1) * vel(ii,jj)
207 15 : write(std_out,*) xcart(ii,jj)
208 : end do
209 : end do
210 :
211 : ! xcart(:,:) = xcart(:,:) + ab_mover%dtion * aa(1) * vel(:,:);
212 :
213 : ! Convert back to xred (reduced coordinates)
214 1 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
215 :
216 : end if ! if (icycle==1)
217 :
218 : !write(std_out,*) 'srkna14 05',jump_end_of_cycle
219 : !##########################################################
220 : !### 05. Compute the next values (Only for extra cycles)
221 :
222 15 : if (icycle>1) then
223 :
224 70 : do ii=1,ab_mover%natom
225 238 : do jj=1,3
226 224 : fcart_m(jj,ii) = fcart(jj,ii)/ab_mover%amass(ii)
227 : end do
228 : end do
229 :
230 14 : if (icycle<16)then
231 :
232 : ! Update of velocities and positions
233 238 : vel(:,:) = vel(:,:) + bb(icycle) * ab_mover%dtion * fcart_m(:,:)
234 : xcart(:,:) = xcart(:,:) +&
235 238 : & aa(icycle) * ab_mover%dtion * vel(:,:)
236 : ! Convert xcart_next to xred_next (reduced coordinates)
237 : ! for scfcv
238 : call xcart2xred(ab_mover%natom, rprimd, xcart,&
239 14 : & xred)
240 :
241 : end if ! (ii<16)
242 :
243 : end if ! if (icycle>1)
244 :
245 : !write(std_out,*) 'srkna14 06',jump_end_of_cycle
246 : !##########################################################
247 : !### 06. Compute the next values (Only for the last cycle)
248 :
249 15 : if(jump_end_of_cycle)then
250 0 : skipcycle=.TRUE.
251 : else
252 15 : skipcycle=.FALSE.
253 : end if
254 :
255 : !write(std_out,*) 'srkna14 07',jump_end_of_cycle
256 : !##########################################################
257 : !### 07. Update the history with the prediction
258 :
259 : !Increase indexes
260 15 : hist%ihist = abihist_findIndex(hist,+1)
261 :
262 : !Fill the history with the variables
263 : !xred, acell, rprimd, vel
264 15 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
265 255 : hist%vel(:,:,hist%ihist)=vel(:,:)
266 15 : ihist_prev = abihist_findIndex(hist,-1)
267 15 : hist%time(hist%ihist)=hist%time(ihist_prev)+ab_mover%dtion
268 :
269 : end subroutine pred_srkna14
270 : !!***
271 :
272 : end module m_pred_srkhna14
273 : !!***
|
