Branch data 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 : : !!***
|
