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