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