Line data Source code
1 : !!****m* ABINIT/m_pred_langevin
2 : !! NAME
3 : !! m_pred_langevin
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_langevin
22 :
23 : use defs_basis
24 : use m_abicore
25 : use m_abimover
26 : use m_abihist
27 :
28 : use m_numeric_tools, only : uniformrandom
29 : use m_geometry, only : xcart2xred, xred2xcart, metric
30 : use m_results_gs , only : results_gs_type
31 :
32 : implicit none
33 :
34 : private
35 : !!***
36 :
37 : public :: pred_langevin
38 : !!***
39 :
40 : contains
41 : !!***
42 :
43 : !!****f* ABINIT/pred_langevin
44 : !! NAME
45 : !! pred_langevin
46 : !!
47 : !! FUNCTION
48 : !! Ionmov predictors (9) Langevin dynamics algorithm
49 : !!
50 : !! IONMOV 9:
51 : !! Uses a Langevin dynamics algorithm :
52 : !! see J. Chelikowsky, J. Phys. D : Appl Phys. 33(2000)R33 [[cite:Chelikowsky2000]]
53 : !!
54 : !! INPUTS
55 : !! ab_mover <type(abimover)> : Datatype with all the information
56 : !! needed by the preditor
57 : !! itime : Index of the present iteration
58 : !! ntime : Maximal number of iterations
59 : !! icycle : Index of the present cycle
60 : !! ncycle : Maximal number of cycles
61 : !! zDEBUG : if true print some debugging information
62 : !!
63 : !! OUTPUT
64 : !!
65 : !! SIDE EFFECTS
66 : !! hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
67 : !!
68 : !! SOURCE
69 :
70 20 : subroutine pred_langevin(ab_mover,hist,icycle,itime,ncycle,ntime,zDEBUG,iexit,skipcycle)
71 :
72 : !Arguments ------------------------------------
73 : !scalars
74 : type(abimover),intent(in) :: ab_mover
75 : type(abihist),intent(inout) :: hist
76 : integer,intent(in) :: itime
77 : integer,intent(in) :: ntime
78 : integer,intent(in) :: iexit
79 : integer,intent(in) :: icycle
80 : integer,intent(inout) :: ncycle
81 : logical,intent(in) :: zDEBUG
82 : logical,intent(out) :: skipcycle
83 :
84 : !Local variables-------------------------------
85 : !scalars
86 : integer :: ii,kk,iatom,idim,iatom1,iatom2,itypat,idum=5,ihist_prev,mcfac
87 : real(dp) :: ucvol,ucvol_next
88 : real(dp),parameter :: v2tol=tol8
89 : real(dp) :: etotal,rescale_vel,ran_num1,ran_num2
90 : real(dp) :: ktemp,dist,distx,disty,distz,maxp1,maxp2,v2nose
91 : real(dp) :: sig_gauss,delxi
92 : logical :: jump_end_of_cycle=.FALSE.
93 : character(len=5000) :: message
94 : !arrays
95 20 : integer,allocatable :: imax_perm(:)
96 : real(dp),allocatable,save :: max_perm(:),pot_perm(:)
97 : real(dp),allocatable,save :: ran_force(:,:),lang_force(:,:)
98 : real(dp),allocatable,save :: fcart_mold(:,:),fcart_m(:,:)
99 :
100 : real(dp) :: acell(3),acell_next(3)
101 : real(dp) :: rprim(3,3),rprimd(3,3),rprimd_next(3,3),rprim_next(3,3)
102 : real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3)
103 40 : real(dp) :: fcart(3,ab_mover%natom)
104 40 : real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
105 40 : real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
106 20 : real(dp) :: vel(3,ab_mover%natom)
107 : real(dp) :: strten(6)
108 :
109 : !***************************************************************************
110 : !Beginning of executable session
111 : !***************************************************************************
112 :
113 20 : if(iexit/=0)then
114 2 : ABI_SFREE(pot_perm)
115 2 : ABI_SFREE(max_perm)
116 : ABI_SFREE(imax_perm)
117 2 : ABI_SFREE(ran_force)
118 2 : ABI_SFREE(lang_force)
119 2 : ABI_SFREE(fcart_mold)
120 2 : ABI_SFREE(fcart_m)
121 : return
122 : end if
123 :
124 18 : jump_end_of_cycle=.FALSE.
125 :
126 : !write(std_out,*) 'langevin 02',jump_end_of_cycle
127 : !##########################################################
128 : !### 02. Allocate the arrays
129 : !### These arrays could be allocated from a previus
130 : !### dataset that exit before itime==ntime
131 :
132 18 : if(itime==1)then
133 2 : ABI_SFREE(pot_perm)
134 2 : ABI_SFREE(max_perm)
135 : ABI_SFREE(imax_perm)
136 2 : ABI_SFREE(ran_force)
137 2 : ABI_SFREE(lang_force)
138 2 : ABI_SFREE(fcart_mold)
139 2 : ABI_SFREE(fcart_m)
140 : end if
141 :
142 18 : if (.not.allocated(pot_perm)) then
143 6 : ABI_MALLOC(pot_perm,(ab_mover%natom))
144 : end if
145 18 : if (.not.allocated(max_perm)) then
146 6 : ABI_MALLOC(max_perm,(ab_mover%ntypat))
147 : end if
148 : if (.not.allocated(imax_perm)) then
149 54 : ABI_MALLOC(imax_perm,(ab_mover%ntypat))
150 : end if
151 18 : if (.not.allocated(ran_force)) then
152 6 : ABI_MALLOC(ran_force,(3,ab_mover%natom))
153 : end if
154 18 : if (.not.allocated(lang_force)) then
155 6 : ABI_MALLOC(lang_force,(3,ab_mover%natom))
156 : end if
157 18 : if (.not.allocated(fcart_mold)) then
158 6 : ABI_MALLOC(fcart_mold,(3,ab_mover%natom))
159 : end if
160 18 : if (.not.allocated(fcart_m)) then
161 6 : ABI_MALLOC(fcart_m,(3,ab_mover%natom))
162 : end if
163 :
164 : !write(std_out,*) 'langevin 03',jump_end_of_cycle
165 : !##########################################################
166 : !### 03. Obtain the present values from the history
167 :
168 18 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
169 :
170 198 : fcart(:,:)=hist%fcart(:,:,hist%ihist)
171 126 : strten(:) =hist%strten(:,hist%ihist)
172 198 : vel(:,:) =hist%vel(:,:,hist%ihist)
173 18 : etotal =hist%etot(hist%ihist)
174 :
175 18 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
176 72 : do ii=1,3
177 234 : rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
178 : end do
179 18 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
180 :
181 : !write(std_out,*) 'langevin 04',jump_end_of_cycle
182 : !##########################################################
183 : !### 04. Compute the next values (Only for the first cycle)
184 :
185 18 : if (icycle==1) then
186 :
187 : ! The temperature is linear between initial and final values
188 : ! It is here converted from Kelvin to Hartree (kb_HaK)
189 18 : ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK
190 : ! write(std_out,*) 'KTEMP=',ktemp
191 : ! write(std_out,*) 'MDITEMP=',ab_mover%mdtemp(1)
192 : ! write(std_out,*) 'MDFTEMP=',ab_mover%mdtemp(2)
193 : ! write(std_out,*) 'DELAYPERM=',ab_mover%delayperm
194 :
195 :
196 : ! %%% LANGEVIN DYNAMICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
197 :
198 : ! 1. Next values are the present ones
199 18 : acell_next(:)=acell(:)
200 18 : ucvol_next=ucvol
201 18 : rprim_next(:,:)=rprim(:,:)
202 18 : rprimd_next(:,:)=rprimd(:,:)
203 :
204 18 : if (zDEBUG) then
205 0 : write (std_out,*) '1. Next values are the present ones'
206 0 : write(std_out,*) 'RPRIMD'
207 0 : do kk=1,3
208 0 : write(std_out,*) rprimd(:,kk)
209 : end do
210 0 : write(std_out,*) 'RPRIM'
211 0 : do kk=1,3
212 0 : write(std_out,*) rprim(:,kk)
213 : end do
214 0 : write(std_out,*) 'ACELL'
215 0 : write(std_out,*) acell(:)
216 : end if
217 :
218 18 : if(itime==1)then
219 :
220 : ! 2. Compute twice the kinetic energy of the system, called v2nose
221 2 : v2nose=0.0_dp
222 7 : do iatom=1,ab_mover%natom
223 22 : do idim=1,3
224 20 : v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
225 : end do
226 : end do
227 2 : if (zDEBUG) then
228 0 : write (std_out,*) '2. Compute twice the kinetic energy of the system, called v2nose'
229 0 : write (std_out,*) 'V2NOSE=',v2nose
230 : end if
231 :
232 : ! 3. If there is no kinetic energy, use random numbers
233 2 : if (v2nose<=v2tol) then
234 2 : v2nose=0.0_dp
235 7 : do iatom=1,ab_mover%natom
236 22 : do idim=1,3
237 : ! uniformrandom returns a uniform random deviate between 0.0 and 1.0
238 : ! if it were always 0 or 1, then the following expression
239 : ! would give the requested temperature
240 15 : vel(idim,iatom)=(1.0_dp-2.0_dp*uniformrandom(idum))*sqrt(ktemp/ab_mover%amass(iatom))
241 : ! Recompute v2nose
242 20 : v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
243 : end do
244 : end do
245 : end if
246 :
247 2 : if (zDEBUG) then
248 0 : write (std_out,*) '3. If there is no kinetic energy, use random numbers'
249 0 : write (std_out,*) 'VEL'
250 0 : do kk=1,ab_mover%natom
251 0 : write (std_out,*) vel(:,kk)
252 : end do
253 0 : write (std_out,*) 'V2NOSE=',v2nose
254 : end if
255 :
256 :
257 : ! Now, rescale the velocities to give the proper temperature
258 2 : rescale_vel=sqrt(3.0_dp*ab_mover%natom*(ab_mover%mdtemp(1))*kb_HaK/v2nose)
259 22 : vel(:,:)=vel(:,:)*rescale_vel
260 : ! Recompute v2nose with the rescaled velocities
261 2 : v2nose=0.0_dp
262 7 : do iatom=1,ab_mover%natom
263 22 : do idim=1,3
264 20 : v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
265 : end do
266 : end do
267 : write(message, '(a)' )&
268 2 : & ' Rescaling or initializing velocities to initial temperature'
269 2 : call wrtout(ab_out,message,'COLL')
270 2 : call wrtout(std_out,message,'COLL')
271 : write(message, '(a,D12.5,a,D12.5)' )&
272 2 : & ' --- Scaling factor : ',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1)
273 2 : call wrtout(ab_out,message,'COLL')
274 2 : call wrtout(std_out,message,'COLL')
275 : write(message, '(a,D12.5)' )&
276 2 : & ' --- Effective temperature',v2nose/3.0_dp/(kb_HaK*ab_mover%natom)
277 2 : call wrtout(ab_out,message,'COLL')
278 2 : call wrtout(std_out,message,'COLL')
279 : ! end if itime==0
280 : end if
281 :
282 : ! This section is devoted to the optional atom permutation (JYR 001114)
283 : ! Two input variables are needed
284 : ! ab_mover%delayperm : is the interval (in time steps) at which
285 : ! atoms are tentatively permuted
286 : ! default value could be 0
287 : ! ab_mover%signperm : is the type of bias for the permutation
288 : ! +1 to favor alternation of species
289 : ! -1 to favor segregation
290 :
291 : ! Force no permutation at initial step
292 18 : if (itime/=1 .and. ab_mover%delayperm/=0 .and. ab_mover%ntypat>2) then
293 0 : if (mod(itime-1,ab_mover%delayperm)==0) then
294 : ! Try commutation of atoms.
295 0 : write(message, '(a)')' Attempt of commutation '
296 0 : call wrtout(ab_out,message,'COLL')
297 0 : call wrtout(std_out,message,'COLL')
298 : ! Compute a 'permutation potential'
299 0 : do iatom=1,ab_mover%natom
300 0 : pot_perm(iatom)=0.0_dp
301 0 : do iatom1=1,ab_mover%natom
302 0 : if (iatom1.ne.iatom) then
303 0 : distx=xcart(1,iatom)-xcart(1,iatom1)
304 0 : distx=distx-acell(1)*nint(distx/acell(1))
305 0 : disty=xcart(2,iatom)-xcart(2,iatom1)
306 0 : disty=disty-acell(2)*nint(disty/acell(2))
307 0 : distz=xcart(3,iatom)-xcart(3,iatom1)
308 0 : distz=distz-acell(3)*nint(distz/acell(3))
309 : ! Here we count each atom below 2 angstr as 1, could be customized
310 0 : dist=sqrt(distx*distx+disty*disty+distz*distz)/3.7807
311 0 : write(std_out,*) iatom,iatom1,dist
312 0 : if (ab_mover%typat(iatom).ne.ab_mover%typat(iatom1)) then
313 : mcfac=-1
314 : else
315 0 : mcfac=1
316 : end if
317 0 : if (dist<1.0_dp) dist=1.0_dp
318 : pot_perm(iatom)=pot_perm(iatom)+mcfac*(ab_mover%signperm)*1.0_dp&
319 0 : & /exp(log(dist)*6.0_dp)
320 : end if
321 : end do
322 : end do
323 0 : write(std_out,*) ' Perm_pot ',pot_perm(:)
324 : ! write(message, '(a,10f12.5)' )' Perm_pot ',&
325 : ! & (pot_perm(iatom1),iatom1=1,ab_mover%natom)
326 : ! call wrtout(ab_out,message,'COLL')
327 : ! call wrtout(std_out,message,'COLL')
328 :
329 : ! Find the two atoms, of different types, with the highest perm_pot
330 0 : max_perm(:)=-1.0d9
331 0 : do iatom=1,ab_mover%natom
332 0 : if (pot_perm(iatom) > max_perm(ab_mover%typat(iatom))) then
333 0 : max_perm(ab_mover%typat(iatom))=pot_perm(iatom)
334 0 : imax_perm(ab_mover%typat(iatom))=iatom
335 : end if
336 : end do
337 :
338 0 : if(zDEBUG)then
339 : ! write(message, '(a,10f12.5)' )' max_Perm ',&
340 : ! & (max_perm(itypat),itypat=1,ab_mover%ntypat)
341 : ! call wrtout(std_out,message,'COLL')
342 : ! write(message, '(a,10i12)' )' imax_Perm ',&
343 : ! & (imax_perm(itypat),itypat=1,ab_mover%ntypat)
344 : ! call wrtout(std_out,message,'COLL')
345 0 : write(std_out,*) 'NTYPAT',ab_mover%ntypat
346 0 : write(message, '(a,10f12.5)' )' max_Perm ',&
347 0 : & (max_perm(:))
348 0 : call wrtout(std_out,message,'COLL')
349 0 : write(message, '(a,10i12)' )' imax_Perm ',&
350 0 : & (imax_perm(:))
351 0 : call wrtout(std_out,message,'COLL')
352 : end if
353 :
354 : ! Loop and keep the 2 largest values
355 0 : if (max_perm(1)>max_perm(2)) then
356 0 : maxp1=max_perm(1)
357 0 : maxp2=max_perm(2)
358 0 : iatom1=imax_perm(1)
359 0 : iatom2=imax_perm(2)
360 : else
361 0 : maxp1=max_perm(2)
362 0 : maxp2=max_perm(1)
363 0 : iatom1=imax_perm(2)
364 0 : iatom2=imax_perm(1)
365 : end if
366 :
367 0 : do itypat=3,ab_mover%ntypat
368 0 : if (max_perm(itypat)>maxp1) then
369 0 : maxp2=maxp1
370 0 : iatom2=iatom1
371 0 : maxp1=max_perm(itypat)
372 0 : iatom1=imax_perm(itypat)
373 0 : else if (max_perm(itypat)>maxp2) then
374 0 : maxp2=max_perm(itypat)
375 0 : iatom2=imax_perm(itypat)
376 : end if
377 : end do
378 0 : write(message, '(2(a,i5))' )' Will commute atom...',iatom1,'...of type ',&
379 0 : & ab_mover%typat(iatom1)
380 0 : call wrtout(ab_out,message,'COLL')
381 0 : call wrtout(std_out,message,'COLL')
382 0 : write(message, '(2(a,i5))' )' with atom...',iatom2,'...of type ',&
383 0 : & ab_mover%typat(iatom2)
384 0 : call wrtout(ab_out,message,'COLL')
385 0 : call wrtout(std_out,message,'COLL')
386 :
387 : ! Commute the atoms positions
388 0 : distx=xcart(1,iatom1)
389 0 : disty=xcart(2,iatom1)
390 0 : distz=xcart(3,iatom1)
391 0 : xcart(1,iatom1)=xcart(1,iatom2)
392 0 : xcart(2,iatom1)=xcart(2,iatom2)
393 0 : xcart(3,iatom1)=xcart(3,iatom2)
394 0 : xcart(1,iatom2)=distx
395 0 : xcart(2,iatom2)=disty
396 0 : xcart(3,iatom2)=distz
397 : ! Convert back to xred (reduced coordinates)
398 0 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
399 :
400 : end if ! if(mod(itime,ab_mover%delayperm)==0)
401 :
402 : else
403 : ! write(std_out,*) "I will jump to the end of cycle"
404 18 : jump_end_of_cycle=.TRUE.
405 :
406 : end if ! if (itime/=1 .and. ab_mover%delayperm/=0 .and. ab_mover%ntypat>2)
407 : ! End of the commutation section
408 :
409 : end if ! if (icycle==1)
410 :
411 18 : ABI_SFREE(imax_perm)
412 :
413 : !write(std_out,*) 'langevin 05',jump_end_of_cycle
414 : !##########################################################
415 : !### 05. Compute the next values (Only for extra cycles)
416 :
417 18 : if (icycle>1) then
418 :
419 : ! write(std_out,*) "Entering internal cycle 2",icycle
420 :
421 : ! If the energy computed is higher than the current
422 : ! (etotal_temp) we have to discard the changes
423 : ! and compute again
424 :
425 : ! write(std_out,*) ch10
426 : ! write(std_out,*) 'EVALUATION FORCES',etotal,hist%etot(abihist_findIndex(hist,-1))
427 : ! write(std_out,*) ch10
428 :
429 : ! This is the worst case (2 evaluations of SCFCV)
430 0 : ihist_prev = abihist_findIndex(hist,-1)
431 0 : if (etotal>hist%etot(ihist_prev).and.icycle==2) then
432 :
433 : ! Discard the changes
434 0 : acell(:) =hist%acell(:,ihist_prev)
435 0 : rprimd(:,:)=hist%rprimd(:,:,ihist_prev)
436 0 : xred(:,:) =hist%xred(:,:,ihist_prev)
437 0 : fcart(:,:) =hist%fcart(:,:,ihist_prev)
438 0 : strten(:) =hist%strten(:,ihist_prev)
439 0 : vel(:,:) =hist%vel(:,:,ihist_prev)
440 0 : etotal =hist%etot(ihist_prev)
441 0 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
442 :
443 : ! distx=xcart(1,iatom1)
444 : ! disty=xcart(2,iatom1)
445 : ! distz=xcart(3,iatom1)
446 : ! xcart(1,iatom1)=xcart(1,iatom2)
447 : ! xcart(2,iatom1)=xcart(2,iatom2)
448 : ! xcart(3,iatom1)=xcart(3,iatom2)
449 : ! xcart(1,iatom2)=distx
450 : ! xcart(2,iatom2)=disty
451 : ! xcart(3,iatom2)=distz
452 :
453 : ! Convert back to xred (reduced coordinates)
454 : ! call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
455 0 : write(message, '(a)' )' Commutation unsuccessful, recomputing the forces'
456 0 : call wrtout(ab_out,message,'COLL')
457 0 : call wrtout(std_out,message,'COLL')
458 :
459 : ! This is the best case (only 1 evaluation of SCFCV)
460 : else
461 :
462 0 : write(message, '(a)')' Commutation successful ! Going on'
463 0 : call wrtout(ab_out,message,'COLL')
464 0 : call wrtout(std_out,message,'COLL')
465 :
466 : ! In thisc case we do not need to compute again SCFCV
467 : ! We avoid the second iteration on ii
468 0 : jump_end_of_cycle=.TRUE.
469 :
470 : end if ! etotal > etotal_temp
471 :
472 : end if ! if (icycle>1)
473 :
474 : !write(std_out,*) 'langevin 06',jump_end_of_cycle
475 : !##########################################################
476 : !### 06. Compute the next values (Only for the last cycle)
477 :
478 18 : if(jump_end_of_cycle)then
479 : ! icycle=ncycle
480 : ! write(std_out,*) 'This is the last cycle, avoid the others and continue'
481 18 : skipcycle=.TRUE.
482 : else
483 0 : skipcycle=.FALSE.
484 : end if
485 :
486 18 : if ((icycle==ncycle).OR.(skipcycle)) then
487 :
488 : ! write(std_out,*) 'ENTERING THE FINAL PART',icycle,ncycle
489 :
490 : ! Specific to Langevin dynamics
491 : ! Initialize an array of random forces
492 : ! No random force at itime=0
493 : ! if (itime==0) then
494 18 : if (itime<0) then
495 :
496 0 : ran_force(:,:)=0.0_dp
497 :
498 : else
499 :
500 63 : do iatom=1,ab_mover%natom
501 : ! sig_gauss is the std deviation of the random distribution
502 45 : sig_gauss=sqrt(2.0_dp*(ab_mover%friction)*ab_mover%amass(iatom)*ktemp)
503 : ! write(std_out,*) 'sig_gauss=',sig_gauss
504 : ! write(std_out,*) 'friction=',ab_mover%friction
505 : ! write(std_out,*) 'ktemp=',ktemp
506 198 : do idim=1,3
507 135 : delxi=2.0_dp
508 306 : do while (delxi >= 1.0_dp)
509 171 : ran_num1=2.0_dp*uniformrandom(idum)-1.0_dp
510 171 : ran_num2=2.0_dp*uniformrandom(idum)-1.0_dp
511 171 : delxi=ran_num1*ran_num1+ran_num2*ran_num2
512 : ! write(std_out,*) delxi,ran_num1,ran_num2
513 : end do
514 : ran_force(idim,iatom)=ran_num1*sqrt(-2.0_dp*log(delxi)/delxi)&
515 180 : & *sig_gauss/sqrt(ab_mover%dtion)
516 : ! write(std_out,*) 'ran_force',ran_force(idim,iatom)
517 : end do
518 : end do
519 :
520 18 : if (zDEBUG) then
521 0 : write (std_out,*) '4. Different forces computed'
522 0 : write (std_out,*) 'RAN_FORCE'
523 0 : do kk=1,ab_mover%natom
524 0 : write (std_out,*) ran_force(:,kk)
525 : end do
526 : end if
527 :
528 :
529 18 : if(zDEBUG)then
530 : ! The distribution should be gaussian
531 0 : delxi=0.0_dp
532 0 : do iatom=1,ab_mover%natom
533 0 : do idim=1,3
534 0 : delxi=delxi+(ran_force(idim,iatom)*ab_mover%dtion)**2
535 : end do
536 : end do
537 0 : delxi=delxi/(3.0_dp*ab_mover%natom)
538 0 : write(message, '(2(a,es22.14))' )' variance =',delxi,' asked =',&
539 0 : & 2.0_dp*(ab_mover%friction)*ab_mover%amass(2)*ktemp*ab_mover%dtion
540 0 : call wrtout(std_out,message,'COLL')
541 : end if
542 : ! end if itime\=0
543 :
544 : end if
545 :
546 : ! zDEBUG
547 : ! write(message, '(a)' )' after initializing ran_force'
548 : ! call wrtout(ab_out,message,'COLL')
549 : ! call wrtout(std_out,message,'COLL')
550 : ! ENDzDEBUG
551 :
552 63 : do iatom=1,ab_mover%natom
553 198 : do idim=1,3
554 135 : fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom)
555 180 : ran_force(idim,iatom)=ran_force(idim,iatom)/ab_mover%amass(iatom)
556 : end do
557 : end do
558 198 : lang_force(:,:)=ran_force(:,:)-(ab_mover%friction)*vel(:,:)+fcart_m(:,:)
559 :
560 18 : if (zDEBUG) then
561 0 : write (std_out,*) '4. Different forces computed'
562 0 : write (std_out,*) 'FCART_M'
563 0 : do kk=1,ab_mover%natom
564 0 : write (std_out,*) fcart_m(:,kk)
565 : end do
566 0 : write (std_out,*) 'RAN_FORCE'
567 0 : do kk=1,ab_mover%natom
568 0 : write (std_out,*) ran_force(:,kk)
569 : end do
570 0 : write (std_out,*) 'LANG_FORCE'
571 0 : do kk=1,ab_mover%natom
572 0 : write (std_out,*) lang_force(:,kk)
573 : end do
574 : end if
575 :
576 : ! zDEBUG
577 : ! write(message, '(a)' )'before verlet'
578 : ! call wrtout(ab_out,message,'COLL')
579 : ! call wrtout(std_out,message,'COLL')
580 : ! ENDzDEBUG
581 :
582 : ! Compute next atomic coordinates using Verlet algorithm
583 :
584 : ! Impose no change of acell, ucvol, rprim, and rprimd
585 18 : acell_next(:)=acell(:)
586 18 : ucvol_next=ucvol
587 18 : rprim_next(:,:)=rprim(:,:)
588 18 : rprimd_next(:,:)=rprimd(:,:)
589 :
590 : ! Convert input xred (reduced coordinates) to xcart (cartesian)
591 18 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
592 : ! Uses the velocity
593 : !
594 : ! If an atom wants to cross the walls, velocity is reversed.
595 : !
596 63 : do iatom=1,ab_mover%natom
597 198 : do idim=1,3
598 : delxi=xcart(idim,iatom)+ab_mover%dtion*vel(idim,iatom)+ &
599 135 : & 0.5_dp*ab_mover%dtion*ab_mover%dtion*lang_force(idim,iatom)
600 135 : if ( (delxi > (rprimd(idim,idim)+(ab_mover%mdwall)) ) .or. &
601 : & (delxi < - (ab_mover%mdwall) ) ) then
602 0 : vel(idim,iatom)=-vel(idim,iatom)
603 : delxi=xcart(idim,iatom)+ab_mover%dtion*vel(idim,iatom)+ &
604 0 : & 0.5_dp*ab_mover%dtion*ab_mover%dtion*lang_force(idim,iatom)
605 : end if
606 180 : xcart_next(idim,iatom)=delxi
607 : end do
608 : end do
609 198 : xcart(:,:)=xcart_next(:,:)
610 18 : if (zDEBUG) then
611 0 : write (std_out,*) '5. If an atom wants to cross the walls, velocity is reversed.'
612 0 : write (std_out,*) 'XCART'
613 0 : do kk=1,ab_mover%natom
614 0 : write (std_out,*) xcart(:,kk)
615 : end do
616 : end if
617 :
618 : ! Convert back to xred_next (reduced coordinates)
619 18 : call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
620 18 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
621 :
622 18 : if (itime==1) then
623 : ! no old forces are available at first step
624 : ! Simple update of the velocity
625 : ! first compute vel_nexthalf for next steps
626 22 : vel(:,:)=vel(:,:)+ab_mover%dtion*lang_force(:,:)
627 : else
628 : ! case itime /= 0 normal verlet integration
629 176 : vel(:,:)=vel(:,:)+0.5_dp*ab_mover%dtion*(fcart_mold(:,:)+lang_force(:,:))
630 : end if
631 18 : if (zDEBUG) then
632 0 : write (std_out,*) '5. Change velocity with verlet'
633 0 : write (std_out,*) 'VEL'
634 0 : do kk=1,ab_mover%natom
635 0 : write (std_out,*) vel(:,kk)
636 : end do
637 : end if
638 :
639 : ! Store 'current force' as 'old force'
640 198 : fcart_mold(:,:)=lang_force(:,:)
641 :
642 : end if ! if (icycle==ncycle)
643 :
644 : !write(std_out,*) 'langevin 07',jump_end_of_cycle
645 : !##########################################################
646 : !### 07. Update the history with the prediction
647 :
648 : !Increase indexes
649 18 : hist%ihist = abihist_findIndex(hist,+1)
650 :
651 18 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
652 198 : hist%vel(:,:,hist%ihist)=vel(:,:)
653 18 : hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
654 :
655 18 : if (ab_mover%delayperm==0 .or. ab_mover%ntypat<=2) ncycle=1
656 18 : if(itime==ntime-1) ncycle=1
657 :
658 : end subroutine pred_langevin
659 : !!***
660 :
661 : end module m_pred_langevin
662 : !!***
|
