Branch data Line data Source code
1 : : !!****m* ABINIT/m_pred_isothermal
2 : : !! NAME
3 : : !! m_pred_isothermal
4 : : !!
5 : : !! FUNCTION
6 : : !! Ionmov predictors (13) Isothermal integrator
7 : : !! This program is decribed in the following paper
8 : : !! Explicit integrators for extended systems dynamics
9 : : !! Glenn J Martyna et al.
10 : : !! Mol. Phys., 1996, Vol. 87, pp. 1117-1157 [[cite:Martyna1996]]
11 : : !!
12 : : !! COPYRIGHT
13 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, JCC, JYR, SE)
14 : : !! This file is distributed under the terms of the
15 : : !! GNU General Public License, see ~abinit/COPYING
16 : : !! or http://www.gnu.org/copyleft/gpl.txt .
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : :
24 : : #include "abi_common.h"
25 : :
26 : : module m_pred_isothermal
27 : :
28 : : use defs_basis
29 : : use m_errors
30 : : use m_abicore
31 : : use m_abimover
32 : : use m_abihist
33 : : use m_linalg_interfaces
34 : :
35 : : use m_numeric_tools, only : uniformrandom
36 : : use m_geometry, only : mkrdim, xcart2xred, xred2xcart, metric
37 : :
38 : : implicit none
39 : :
40 : : private
41 : : !!***
42 : :
43 : : public :: pred_isothermal
44 : : !!***
45 : :
46 : : contains
47 : : !!***
48 : :
49 : : !!****f* ABINIT/pred_isothermal
50 : : !! NAME
51 : : !! pred_isothermal
52 : : !!
53 : : !! FUNCTION
54 : : !! Ionmov predictors (13) Isothermal integrator
55 : : !!
56 : : !! IONMOV 13:
57 : : !! Reversible integrator of Martyna at al.
58 : : !! The equation of motion of the ions in contact with a thermostat
59 : : !! and a barostat are solved with the algorithm proposed by Martyna,
60 : : !! Tuckermann Tobias and Klein, Mol. Phys., 1996, p. 1117. [[cite:Martyna1996]]
61 : : !! Related parameters : the time step (dtion),
62 : : !! the initial temperature mdtemp(1), the final temperature mdtemp(2),
63 : : !! the number of thermostats (nnos), and the masses of thermostats (qmass).
64 : : !! If optcell=1 or 2, the mass of the barostat (bmass) must be given in addition.
65 : : !!
66 : : !! There are three sub cases according to the value of optcell
67 : : !! optcell=0: isothermal
68 : : !! optcell=1: homogeneous cell fluctuations
69 : : !! optcell=2: full cell fluctuation in addition to temperature control.
70 : : !!
71 : : !! INPUTS
72 : : !! ab_mover <type(abimover)> : Datatype with all the information
73 : : !! needed by the preditor
74 : : !! itime : Index of the present iteration
75 : : !! ntime : Maximal number of iterations
76 : : !! zDEBUG : if true print some debugging information
77 : : !!
78 : : !! SIDE EFFECTS
79 : : !! hist <type(abihist)> : History of positions,forces
80 : : !! acell, rprimd, stresses
81 : : !!
82 : : !! SOURCE
83 : :
84 : 9093 : subroutine pred_isothermal(ab_mover,hist,itime,mttk_vars,ntime,zDEBUG,iexit)
85 : :
86 : : !Arguments ------------------------------------
87 : : !scalars
88 : : type(abimover),intent(in) :: ab_mover
89 : : type(abihist),intent(inout) :: hist
90 : : type(mttk_type),intent(inout) :: mttk_vars
91 : : integer,intent(in) :: itime
92 : : integer,intent(in) :: ntime
93 : : integer,intent(in) :: iexit
94 : : logical,intent(in) :: zDEBUG
95 : :
96 : : !Local variables-------------------------------
97 : : !scalars
98 : : integer :: ii,kk,iatom,idim,idum=5,ierr
99 : : integer,parameter :: lwork=8
100 : : real(dp) :: ucvol,ucvol0,ucvol_next,mttk_aloc,mttk_aloc2,mttk_bloc,ekin
101 : : real(dp) :: massvol=0
102 : : real(dp),parameter :: esh2=one/six,esh4=esh2/20._dp,esh6=esh4/42._dp
103 : : real(dp),parameter :: esh8=esh6/72._dp,nosetol=tol10,v2tol=tol8
104 : : real(dp) :: etotal,rescale_vel,polysh,s1,s2,sigma2,v2gauss,vtest
105 : : real(dp),save :: ktemp,vlogv
106 : : character(len=5000) :: message
107 : : !arrays
108 : : real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:)
109 : :
110 : : real(dp) :: mttk_alc(3),mttk_alc2(3),mttk_blc(3),mttk_psh(3)
111 : : real(dp) :: mttk_tv(3,3),mttk_vt(3,3),mttk_ubox(3,3)
112 : : real(dp) :: mttk_uu(3),mttk_uv(3),mttk_veig(3)
113 : : real(dp) :: acell(3),acell0(3),acell_next(3)
114 : : real(dp) :: rprimd(3,3),rprimd0(3,3),rprim(3,3),rprimd_next(3,3),rprim_next(3,3)
115 : : real(dp) :: gprimd(3,3)
116 : : real(dp) :: gmet(3,3)
117 : : real(dp) :: rmet(3,3)
118 : 18186 : real(dp) :: fcart(3,ab_mover%natom)
119 : 18186 : real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
120 : 18186 : real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
121 : 9093 : real(dp) :: vel(3,ab_mover%natom)
122 : : real(dp) :: strten(6),work(lwork)
123 : :
124 : : !***************************************************************************
125 : : !Beginning of executable session
126 : : !***************************************************************************
127 : :
128 [ + + ]: 9093 : if(iexit/=0)then
129 [ + - ]: 31 : ABI_SFREE(fcart_m)
130 [ + - ]: 31 : ABI_SFREE(vel_nexthalf)
131 : : return
132 : : end if
133 : :
134 : : !write(std_out,*) 'isothermal 01'
135 : : !##########################################################
136 : : !### 01. Debugging and Verbose
137 : :
138 [ - + ]: 9062 : if(zDEBUG)then
139 [ # # ]: 0 : write(std_out,'(a,3a,41a,36a)') ch10,('-',kk=1,3),&
140 [ # # # # : 0 : & 'Debugging and Verbose for pred_isothermal',('-',kk=1,36)
# # ]
141 : 0 : write(std_out,*) 'ionmov: ',13
142 : 0 : write(std_out,*) 'itime: ',itime
143 : : end if
144 : :
145 : : !write(std_out,*) 'isothermal 01'
146 : : !##########################################################
147 : : !### 01. Allocate the vectors vin, vout and hessian matrix
148 : : !### These arrays could be allocated from a previus
149 : : !### dataset that exit before itime==ntime
150 : :
151 [ + + ]: 9062 : if(itime==1)then
152 [ + + ]: 233 : ABI_SFREE(fcart_m)
153 [ + + ]: 233 : ABI_SFREE(vel_nexthalf)
154 : : end if
155 : :
156 [ + + ]: 9062 : if (.not.allocated(fcart_m)) then
157 [ - + + - : 699 : ABI_MALLOC(fcart_m,(3,ab_mover%natom))
+ - + - -
+ - + -
+ ]
158 : : end if
159 [ + + ]: 9062 : if (.not.allocated(vel_nexthalf)) then
160 [ - + + - : 699 : ABI_MALLOC(vel_nexthalf,(3,ab_mover%natom))
+ - + - -
+ - + -
+ ]
161 : : end if
162 : :
163 : : !write(std_out,*) 'isothermal 02'
164 : : !##########################################################
165 : : !### 02. Obtain the present values from the history
166 : :
167 : 9062 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
168 : :
169 [ + + + + ]: 35378022 : fcart(:,:)=hist%fcart(:,:,hist%ihist)
170 [ + + ]: 63434 : strten(:) =hist%strten(:,hist%ihist)
171 [ + + + + ]: 35378022 : vel(:,:) =hist%vel(:,:,hist%ihist)
172 : 9062 : etotal =hist%etot(hist%ihist)
173 : :
174 [ + + ]: 36248 : do ii=1,3
175 [ + + ]: 117806 : rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
176 : : end do
177 : 9062 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
178 : :
179 : 9062 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
180 : :
181 [ - + ]: 9062 : if(zDEBUG)then
182 : 0 : write (std_out,*) 'fcart:'
183 [ # # ]: 0 : do kk=1,ab_mover%natom
184 : 0 : write (std_out,*) fcart(:,kk)
185 : : end do
186 : 0 : write (std_out,*) 'vel:'
187 [ # # ]: 0 : do kk=1,ab_mover%natom
188 : 0 : write (std_out,*) vel(:,kk)
189 : : end do
190 : 0 : write (std_out,*) 'strten:'
191 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
192 : 0 : write (std_out,*) 'etotal:'
193 : 0 : write (std_out,*) etotal
194 : : end if
195 : :
196 : : !Save initial values
197 : 9062 : acell0(:)=acell(:)
198 : 9062 : rprimd0(:,:)=rprimd(:,:)
199 : 9062 : ucvol0=ucvol
200 : :
201 : : !write(std_out,*) 'isothermal 03'
202 : : !##########################################################
203 : : !### 05. Seconde half velocity step
204 : :
205 [ + + ]: 9062 : if (itime>1) then
206 : :
207 : : ! Next Half velocity step
208 [ + + ]: 35316 : do idim=1,3
209 [ + + ]: 25904088 : fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:)
210 : : end do
211 [ + + + + ]: 34500525 : vel(:,:)=vel_nexthalf(:,:)+ab_mover%dtion/two*fcart_m(:,:)
212 : :
213 [ + + ]: 8829 : if (ab_mover%optcell==0) then
214 : : ! Update Thermostat variables and velocity
215 : : call isotemp(ab_mover%amass,ab_mover%dtion,ekin,ab_mover%iatfix,&
216 [ + - - - : 7861 : & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,vel)
- - - + ]
217 [ + + ]: 968 : else if (ab_mover%optcell==1) then
218 : : ! Update Thermostat variables and velocity
219 : : call isopress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,&
220 : : & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,&
221 [ + - - - : 14 : & strten,ab_mover%strtarget,ucvol,vel,vlogv)
- - + - -
- - - - +
- + ]
222 [ + - ]: 954 : else if (ab_mover%optcell==2) then
223 : : ! Next half step for extended variables
224 : : call isostress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,&
225 : : & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,&
226 [ + - - - : 954 : & ab_mover%qmass,strten,ab_mover%strtarget,ucvol,vel)
- - + - -
- - - - +
- + ]
227 : : end if
228 : :
229 [ + + ]: 8829 : if(itime==2) massvol=ekin+etotal
230 : :
231 [ + + ]: 8829 : if (ab_mover%optcell==2) then
232 : : ! Evolution of cell and volume
233 : 954 : acell_next(:)=acell(:)
234 : 954 : ucvol_next=ucvol
235 : : end if
236 : :
237 : 8829 : call mkrdim(acell,rprim,rprimd)
238 : 8829 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
239 : :
240 [ - + ]: 8829 : if(zDEBUG)then
241 : 0 : write(std_out,*) 'Second half velocity step'
242 : 0 : write(std_out,*) 'Cell parameters:'
243 : 0 : write(std_out,*) 'rprimd:'
244 [ # # ]: 0 : do kk=1,3
245 : 0 : write(std_out,*) rprimd(:,kk)
246 : : end do
247 : 0 : write(std_out,*) 'rprim:'
248 [ # # ]: 0 : do kk=1,3
249 : 0 : write(std_out,*) rprim(:,kk)
250 : : end do
251 : 0 : write(std_out,*) 'acell:'
252 : 0 : write(std_out,*) acell(:)
253 : 0 : write(std_out,*) 'Conserved energy:',(ekin+etotal)-massvol,ekin,etotal
254 : 0 : write(std_out,*) 'Volume of unitqry cell (ucvol):',ucvol
255 : : end if
256 : :
257 : : end if ! if (itime>1)
258 : :
259 : : !write(std_out,*) 'isothermal 04'
260 : : !##########################################################
261 : : !### 03. Compute the next values
262 : :
263 : : !The temperature is linear between initial and final values
264 : : !It is here converted from Kelvin to Hartree (kb_HaK)
265 : 9062 : ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK
266 : :
267 [ - + ]: 9062 : if(zDEBUG)then
268 : 0 : write(std_out,*) 'Temperature in Kelvin (ktemp):',ktemp
269 : 0 : write(std_out,*) 'Initial temp (mdtemp(1)):',ab_mover%mdtemp(1)
270 : 0 : write(std_out,*) 'Final temp (mdtemp(2)):',ab_mover%mdtemp(2)
271 : 0 : write(std_out,*) 'Delay for atom permutation (delayperm)',ab_mover%delayperm
272 : 0 : write(std_out,*) 'dtion',ab_mover%dtion
273 : 0 : write(std_out,*) 'nnos:', ab_mover%nnos
274 : 0 : write(std_out,*) 'qmass', ab_mover%qmass(:)
275 : 0 : write(std_out,*) 'bmass',ab_mover%bmass
276 : : end if
277 : :
278 [ + + ]: 9062 : if(itime==1) then
279 [ + + + + : 962 : mttk_vars%glogs(:)=zero; mttk_vars%vlogs(:)=zero; mttk_vars%xlogs(:)=zero
+ + ]
280 [ + + + + ]: 3029 : mttk_vars%vboxg(:,:)=zero
281 : 233 : vlogv=zero
282 : : ! v2gauss is twice the kinetic energy
283 : 233 : v2gauss=0.0_dp
284 [ + + ]: 219549 : do iatom=1,ab_mover%natom
285 [ + + ]: 877497 : do idim=1,3
286 : 877264 : v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
287 : : end do
288 : : end do
289 : :
290 : : ! If there is no kinetic energy
291 [ + - ]: 233 : if (v2gauss<=v2tol.and.itime==1) then
292 : : ! Maxwell-Boltzman distribution
293 : 219549 : v2gauss=zero
294 : 219549 : vtest=zero
295 [ + + ]: 219549 : do iatom=1,ab_mover%natom
296 [ + + ]: 877497 : do idim=1,3
297 : 657948 : vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum))
298 : 877264 : vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum)))
299 : : end do
300 : : end do
301 : :
302 : : ! Get rid of center-of-mass velocity
303 [ + + ]: 219549 : s1=sum(ab_mover%amass(:))
304 [ + + ]: 932 : do idim=1,3
305 [ + + ]: 658647 : s2=sum(ab_mover%amass(:)*vel(idim,:))
306 [ + + ]: 658880 : vel(idim,:)=vel(idim,:)-s2/s1
307 : : end do
308 : :
309 : : ! Recompute v2gauss
310 [ + + ]: 219549 : do iatom=1,ab_mover%natom
311 [ + + ]: 877497 : do idim=1,3
312 : 657948 : v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
313 : 877264 : vtest=vtest+vel(idim,iatom)/(3._dp*ab_mover%natom)
314 : : end do
315 : : end do
316 : :
317 : : ! Now rescale the velocities to give the exact temperature
318 : 233 : rescale_vel=sqrt(3._dp*ab_mover%natom*kb_HaK*ab_mover%mdtemp(1)/v2gauss)
319 [ + + + + ]: 877497 : vel(:,:)=vel(:,:)*rescale_vel
320 : :
321 : : ! Recompute v2gauss with the rescaled velocities
322 : : v2gauss=zero
323 [ + + ]: 219549 : do iatom=1,ab_mover%natom
324 [ + + ]: 877497 : do idim=1,3
325 : 877264 : v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
326 : : end do
327 : : end do
328 : :
329 : : ! Compute the variance and print
330 : 233 : sigma2=(v2gauss/(3._dp*ab_mover%natom)-ab_mover%amass(1)*vtest**2)/kb_HaK
331 : :
332 [ - + ]: 233 : if (zDEBUG)then
333 : : write(message, '(a)' )&
334 : 0 : & ' Rescaling or initializing velocities to initial temperature'
335 : 0 : call wrtout(std_out,message,'COLL')
336 : : write(message, '(a,d12.5,a,D12.5)' )&
337 : 0 : & ' --- Scaling factor :',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1)
338 : 0 : call wrtout(std_out,message,'COLL')
339 : : write(message, '(a,d12.5,a,D12.5)' )&
340 : 0 : & ' --- Effective temperature',v2gauss/(3*ab_mover%natom*kb_HaK),' From variance', sigma2
341 : 0 : call wrtout(std_out,message,'COLL')
342 : : end if
343 : :
344 : : end if !(v2gauss<=v2tol.and.itime==1)
345 : : end if !(itime==1)
346 : :
347 : : !XG070613 : Do not take away the following line , seems needed for the pathscale compiler
348 : :
349 [ - + ]: 9062 : if (zDEBUG) write(std_out,*) 'vboxg',mttk_vars%vboxg(:,:)
350 : :
351 : :
352 : : !write(std_out,*) 'isothermal 05'
353 : : !##########################################################
354 : : !### 03. First half velocity step
355 : :
356 : : !write(std_out,*) 'FIRST HALF VELOCITY STEP',ucvol
357 : : !write(std_out,*) 'OPTCELL option selected:',ab_mover%optcell
358 : : !write(std_out,*) 'RPRIMD'
359 : : !do kk=1,3
360 : : !write(std_out,*) rprimd(:,kk)
361 : : !end do
362 : : !write(std_out,*) 'RPRIM'
363 : : !do kk=1,3
364 : : !write(std_out,*) rprim(:,kk)
365 : : !end do
366 : : !write(std_out,*) 'ACELL'
367 : : !write(std_out,*) acell(:)
368 : :
369 : : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
370 : : !%%% BEGIN sub case optcell=0 Isothermal Ensemble
371 : : !%%%
372 [ + + ]: 9062 : if(ab_mover%optcell==0) then
373 : : ! There is no evolution of cell
374 : 8068 : acell_next(:)=acell(:)
375 : 8068 : ucvol_next=ucvol
376 : 8068 : rprim_next(:,:)=rprim(:,:)
377 : 8068 : rprimd_next(:,:)=rprimd(:,:)
378 : : ! Update Thermostat variables and scale velocitie
379 : : call isotemp(ab_mover%amass,ab_mover%dtion,ekin,ab_mover%iatfix,&
380 [ + - - - : 8068 : & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,vel)
- - - + ]
381 : :
382 : : ! Half velocity step
383 [ + + ]: 32272 : do idim=1,3
384 [ + + ]: 25966768 : fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:)
385 : : end do
386 [ + + + + ]: 34587396 : vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:)
387 : : ! New positions
388 : : ! Convert input xred (reduced coordinates) to xcart (cartesian)
389 : 8068 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
390 [ + + + + ]: 34587396 : xcart_next(:,:)=xcart(:,:)+vel_nexthalf(:,:)*ab_mover%dtion
391 : : ! Convert back to xred (reduced coordinates)
392 : 8068 : call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
393 : : ! %%%
394 : : ! %%% END sub case optcell=0 Isothermal Ensemble
395 : : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
396 : :
397 : : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
398 : : ! %%% BEGIN sub case optcell=1 Isothermal-Isenthalpic
399 : : ! %%% Ensemble (homogeneous cell deformation)
400 : : ! %%%
401 [ + + ]: 994 : else if (ab_mover%optcell==1) then
402 : : ! Only homogeneous evolution of cell
403 : : ! Evolution of cell we keep rprim constant
404 : 16 : rprim_next(:,:)=rprim(:,:)
405 : : ! Update Thermostat variables and velocity
406 : : call isopress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,&
407 : : & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,&
408 [ + - - - : 16 : & strten,ab_mover%strtarget,ucvol,vel,vlogv)
- - + - -
- - - - +
- + ]
409 : :
410 : : ! Half velocity step
411 [ + + ]: 64 : do idim=1,3
412 [ + + ]: 256 : fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:)
413 : : end do
414 [ + + + + ]: 272 : vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:)
415 : : ! New positions
416 : 16 : mttk_aloc=exp(ab_mover%dtion/two*vlogv)
417 : 16 : mttk_aloc2=(vlogv*ab_mover%dtion/two)**2
418 : 16 : polysh=(((esh8*mttk_aloc2+esh6)*mttk_aloc2+esh4)*mttk_aloc2+esh2)*mttk_aloc2+one
419 : 16 : mttk_bloc=mttk_aloc*polysh*ab_mover%dtion
420 : : ! Convert input xred (reduced coordinates) to xcart (cartesian)
421 : 16 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
422 [ + + + + ]: 272 : xcart_next(:,:)=xcart(:,:)*mttk_aloc**2+vel_nexthalf(:,:)*mttk_bloc
423 : : ! Update the volume and related quantities
424 [ + + ]: 64 : acell_next(:)=acell(:)*exp(ab_mover%dtion*vlogv)
425 : : ! ucvol=ucvol*exp(ab_mover%dtion*vlogv)
426 : 16 : call mkrdim(acell_next,rprim,rprimd_next)
427 : 16 : call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol_next)
428 : : ! Convert back to xred (reduced coordinates)
429 : 16 : call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next)
430 : : ! Computation of the forces for the new positions
431 : : ! Compute DFT forces (big loop)
432 : :
433 : : ! COMMENTED
434 : : ! This should be in mover.F90
435 : :
436 : : ! ! If metric has changed since the initialization, update the Ylm's
437 : : ! if (ab_mover%optcell/=0.and.psps%useylm==1.and.itime>1)then
438 : : ! option=0;if (ab_mover%iscf>0) option=1
439 : : ! call initylmg(gprimd,kg,ab_mover%kptns,ab_mover%mkmem,mpi_enreg,psps%mpsang,ab_mover%mpw,ab_mover%nband,ab_mover%nkpt,&
440 : : ! & npwarr,ab_mover%nsppol,option,rprimd_next,ylm,ylmgr)
441 : : ! end if
442 : :
443 : :
444 : : ! %%%
445 : : ! %%% END sub case optcell=1 Isothermal-Isenthalpic
446 : : ! %%% Ensemble (homogeneous cell deformation)
447 : : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
448 : :
449 : : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
450 : : ! %%% BEGIN sub case optcell=2 Isothermal-Isenthalpic
451 : : ! %%% Ensemble (full cell deformation)
452 : : ! %%%
453 [ + - ]: 978 : else if (ab_mover%optcell==2) then
454 : 978 : acell_next=acell
455 : : ! Fisrt half step for extended variables
456 : : call isostress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,&
457 : : & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,&
458 [ + - - - : 978 : & ab_mover%qmass,strten,ab_mover%strtarget,ucvol,vel)
- - + - -
- - - - +
- + ]
459 : : ! Half velocity step
460 [ + + ]: 3912 : do idim=1,3
461 [ + + ]: 595944 : fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:)
462 : : end do
463 [ + + + + ]: 790354 : vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:)
464 : : ! Convert input xred (reduced coordinates) to xcart (cartesian)
465 : 978 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
466 : : ! New positions
467 [ + + + + ]: 12714 : mttk_vt(:,:)=mttk_vars%vboxg(:,:)
468 : 978 : call dsyev('V','U',3,mttk_vt,3,mttk_veig,work,lwork,ierr)
469 [ + + + + ]: 12714 : mttk_tv(:,:)=transpose(mttk_vt)
470 [ + + ]: 3912 : mttk_alc(:)=exp(ab_mover%dtion/two*mttk_veig(:))
471 [ + + ]: 3912 : mttk_alc2(:)=(mttk_veig(:)*ab_mover%dtion/two)**2
472 [ + + ]: 3912 : mttk_psh(:)=(((esh8*mttk_alc2(:)+esh6)*mttk_alc2(:)+esh4)*mttk_alc2(:)+esh2)*mttk_alc2(:)+one
473 [ + + ]: 3912 : mttk_blc(:)=mttk_alc(:)*mttk_psh(:)*ab_mover%dtion
474 : : ! Update the positions
475 [ + + ]: 198322 : do iatom=1,ab_mover%natom
476 [ + + + + ]: 2565472 : mttk_uu(:)=matmul(mttk_tv,xcart(:,iatom))
477 [ + + + + : 2565472 : mttk_uv(:)=matmul(mttk_tv,vel_nexthalf(:,iatom))
+ - ]
478 [ + + ]: 789376 : mttk_uu(:)=mttk_uu(:)*mttk_alc(:)**2+mttk_uv(:)*mttk_blc(:)
479 [ + + + + : 3158482 : xcart_next(:,iatom)=matmul(mttk_vt,mttk_uu)
+ + ]
480 : : end do
481 : : ! Update the box (rprimd and rprim)
482 [ + + + + : 39120 : mttk_ubox(:,:)=matmul(mttk_tv,rprimd)
+ + ]
483 [ + + ]: 3912 : do idim=1,3
484 [ + + ]: 12714 : mttk_ubox(:,idim)=mttk_ubox(:,idim)*mttk_alc(:)**2
485 : : end do
486 [ + + + + : 39120 : rprimd_next(:,:)=matmul(mttk_vt,mttk_ubox)
+ + ]
487 [ + + ]: 3912 : do idim=1,3
488 [ + + ]: 12714 : rprim_next(idim,:)=rprimd_next(idim,:)/acell(:)
489 : : end do
490 : : ! Update the volume
491 : 978 : call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol)
492 : : ! Convert back to xred (reduced coordinates)
493 : 978 : call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next)
494 : : ! Computation of the forces for the new positions
495 : :
496 : : ! COMMENTED
497 : : ! This should be in mover.F90
498 : :
499 : : ! ! If metric has changed since the initialization, update the Ylm's
500 : : ! if (ab_mover%optcell/=0.and.psps%useylm==1.and.itime>1)then
501 : : ! option=0;if (ab_mover%iscf>0) option=1
502 : : ! call initylmg(gprimd,kg,ab_mover%kptns,ab_mover%mkmem,mpi_enreg,psps%mpsang,ab_mover%mpw,ab_mover%nband,ab_mover%nkpt,&
503 : : ! & npwarr,ab_mover%nsppol,option,rprimd_next,ylm,ylmgr)
504 : : ! end if
505 : :
506 : : ! %%%
507 : : ! %%% END sub case optcell=2 Isothermal-Isenthalpic
508 : : ! %%% Ensemble (full cell deformation)
509 : : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
510 : : else
511 : : write(message, '(a,i12,a,a)' )&
512 : 0 : & ' Disallowed value for optcell=',ab_mover%optcell,ch10,&
513 : 0 : & ' Allowed values with ionmov==13 : 0 to 2.'
514 : 0 : ABI_BUG(message)
515 : : end if
516 : :
517 : : !write(std_out,*) 'OLD PARAMETERS'
518 : : !write(std_out,*) 'RPRIMD'
519 : : !do kk=1,3
520 : : !write(std_out,*) rprimd(:,kk)
521 : : !end do
522 : : !write(std_out,*) 'RPRIM'
523 : : !do kk=1,3
524 : : !write(std_out,*) rprim(:,kk)
525 : : !end do
526 : : !write(std_out,*) 'ACELL'
527 : : !write(std_out,*) acell(:)
528 : :
529 : : !write(std_out,*) 'NEXT PARAMETERS'
530 : : !write(std_out,*) 'RPRIMD'
531 : : !do kk=1,3
532 : : !write(std_out,*) rprimd_next(:,kk)
533 : : !end do
534 : : !write(std_out,*) 'RPRIM'
535 : : !do kk=1,3
536 : : !write(std_out,*) rprim_next(:,kk)
537 : : !end do
538 : : !write(std_out,*) 'ACELL'
539 : : !write(std_out,*) acell_next(:)
540 : :
541 : :
542 : : !Those are the values store into the history
543 : 9062 : rprim=rprim_next
544 : 9062 : rprimd=rprimd_next
545 [ + + + + ]: 35378022 : xred=xred_next
546 [ + + + + ]: 35378022 : xcart=xcart_next
547 : 9062 : acell=acell_next
548 : :
549 : : !write(std_out,*) 'isothermal 06'
550 : : !##########################################################
551 : : !### 06. Update the history with the prediction
552 : :
553 : : !Increase indexes
554 : 9062 : hist%ihist = abihist_findIndex(hist,+1)
555 : :
556 : : !Fill the history with the variables
557 : : !xred, acell, rprimd, vel
558 : 9062 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
559 [ + + + + ]: 35378022 : hist%vel(:,:,hist%ihist)=vel(:,:)
560 : 9062 : hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
561 : :
562 [ - + ]: 9062 : if(zDEBUG)then
563 : 0 : write (std_out,*) 'fcart:'
564 [ # # ]: 0 : do kk=1,ab_mover%natom
565 : 0 : write (std_out,*) fcart(:,kk)
566 : : end do
567 : 0 : write (std_out,*) 'vel:'
568 [ # # ]: 0 : do kk=1,ab_mover%natom
569 : 0 : write (std_out,*) vel(:,kk)
570 : : end do
571 : 0 : write (std_out,*) 'strten:'
572 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
573 : 0 : write (std_out,*) 'etotal:'
574 : 0 : write (std_out,*) etotal
575 : : end if
576 : :
577 : : end subroutine pred_isothermal
578 : : !!***
579 : :
580 : : !!****f* ABINIT/isotemp
581 : : !! NAME
582 : : !! isotemp
583 : : !!
584 : : !! FUNCTION
585 : : !! performs one half step on isotemp parameters according to Martyna et al.
586 : : !!
587 : : !! INPUTS
588 : : !! amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...)
589 : : !! dtion=
590 : : !! isotemp_data
591 : : !! ktemp
592 : : !! vel
593 : : !!
594 : : !! OUTPUT
595 : : !! Only updates variables
596 : : !!
597 : : !! SIDE EFFECTS
598 : : !! isotemp_data: updates the thermostat parameters
599 : : !! vel=update the velocities
600 : : !!
601 : : !! SOURCE
602 : :
603 [ - + - + ]: 15929 : subroutine isotemp(amass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,qmass,vel)
604 : :
605 : : !Arguments ------------------------------------
606 : : !scalars
607 : : integer,intent(in) :: natom
608 : : integer,intent(in) :: nnos
609 : : real(dp),intent(in) :: dtion,ktemp
610 : : real(dp),intent(out) :: ekin
611 : : type(mttk_type) :: mttk_vars
612 : : !arrays
613 : : real(dp),intent(in) :: amass(natom)
614 : : real(dp),intent(inout) :: vel(3,natom)
615 : : real(dp),intent(in) :: qmass(:)
616 : : integer,intent(in) :: iatfix(:,:)
617 : :
618 : : !Local variables ------------------------------
619 : : !scalars
620 : : integer :: iatom,idir,inos
621 : : real(dp) :: alocal,gnkt,nfree,scale
622 : : !character(len=500) :: message
623 : : !arrays
624 : 15929 : real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:)
625 : :
626 : : !***************************************************************************
627 : : !Beginning of executable session
628 : : !***************************************************************************
629 : :
630 [ - + + - : 47787 : ABI_MALLOC(glogs,(nnos))
+ - + - -
+ - + ]
631 [ - + + - : 31858 : ABI_MALLOC(vlogs,(nnos))
- + - + ]
632 [ - + + - : 31858 : ABI_MALLOC(xlogs,(nnos))
- + - + ]
633 [ + + ]: 31888 : glogs(:)=mttk_vars%glogs(:)
634 [ + + ]: 31888 : vlogs(:)=mttk_vars%vlogs(:)
635 [ + + ]: 31888 : xlogs(:)=mttk_vars%xlogs(:)
636 : 15929 : scale=one
637 : : !Compute the ionic kinetic energy
638 : 15929 : nfree=zero
639 : 15929 : ekin=zero
640 [ + + ]: 17089109 : do iatom=1,natom
641 [ + + ]: 68308649 : do idir=1,3
642 : : ! Warning : the fixing of atomis is implemented in reduced
643 : : ! coordinates, so that this expression is wrong
644 [ + - ]: 68292720 : if (iatfix(idir,iatom) == 0) then
645 : 51219540 : ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2
646 : : ! Counts the degrees of freedom
647 : 51219540 : nfree=nfree+one
648 : : end if
649 : : end do
650 : : end do
651 : 15929 : gnkt=nfree*ktemp
652 : : !Update the forces
653 : 15929 : glogs(1)=(two*ekin-gnkt)/qmass(1)
654 : 15929 : vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
655 [ + + ]: 15959 : do inos=1,nnos-1
656 : 30 : alocal=exp(-dtion/eight*vlogs(nnos+1-inos))
657 : : vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+&
658 : 15959 : & dtion/four*glogs(nnos-inos)*alocal
659 : : end do
660 : : !Update the particle velocities
661 : 15929 : alocal=exp(-dtion/two*vlogs(1))
662 : 15929 : scale=scale*alocal
663 : : !Update the forces
664 : 15929 : glogs(1)=(scale*scale*two*ekin-gnkt)/qmass(1)
665 : : !Update the thermostat positions
666 [ + + ]: 31888 : do inos=1,nnos
667 : 31888 : xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two
668 : : end do
669 : : !Update the thermostat velocities
670 [ + + ]: 15959 : do inos=1,nnos-1
671 : 30 : alocal=exp(-dtion/eight*vlogs(inos+1))
672 : 30 : vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal
673 : 15959 : glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1)
674 : : end do
675 : 15929 : vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
676 [ + + + + ]: 68308649 : vel(:,:)=vel(:,:)*scale
677 : : !Compute the ionic kinetic energy
678 : 15929 : ekin=zero
679 [ + + ]: 17089109 : do iatom=1,natom
680 [ + + ]: 68308649 : do idir=1,3
681 : : ! Warning : the fixing of atomis is implemented in reduced
682 : : ! coordinates, so that this expression is wrong
683 [ + - ]: 68292720 : if (iatfix(idir,iatom) == 0) then
684 : 51219540 : ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2
685 : : end if
686 : : end do
687 : : end do
688 : : !Compute the thermostat kinetic energy and add it to the ionic one
689 : 15929 : ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*nfree*ktemp
690 [ + + ]: 15959 : do inos=2,nnos
691 : 15959 : ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp
692 : : end do
693 [ + + ]: 31888 : mttk_vars%glogs(:)=glogs(:)
694 [ + + ]: 31888 : mttk_vars%vlogs(:)=vlogs(:)
695 [ + + ]: 31888 : mttk_vars%xlogs(:)=xlogs(:)
696 : 15929 : ABI_FREE(glogs)
697 : 15929 : ABI_FREE(vlogs)
698 : 15929 : ABI_FREE(xlogs)
699 : : !DEBUG
700 : : !write(std_out,*)'ekin added',half*qmass(1)*vlogs(1)**2,xlogs(1)*(nfree)*ktemp
701 : : !ENDEBUG
702 : :
703 : 15929 : end subroutine isotemp
704 : : !!***
705 : :
706 : : !----------------------------------------------------------------------
707 : :
708 : : !!****f* ABINIT/isopress
709 : : !! NAME
710 : : !! isopress
711 : : !!
712 : : !! FUNCTION
713 : : !! performs one half step on isopress parameters according to Martyna et al.
714 : : !!
715 : : !! INPUTS
716 : : !! amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...)
717 : : !! dtion= ionic time step
718 : : !! isotemp_data
719 : : !! ktemp
720 : : !! press= current pressure of the system
721 : : !! prtarget= target pressure
722 : : !! ucvol= unit cell volume
723 : : !! vel= current velocity
724 : : !!
725 : : !! OUTPUT
726 : : !! Only updates variables
727 : : !!
728 : : !! SIDE EFFECTS
729 : : !! isotemp_data: updates the thermostat parameters (saved variables: bouh !)
730 : : !! vel=update the velocities
731 : : !!
732 : : !! SOURCE
733 : :
734 [ - + ]: 30 : subroutine isopress(amass,bmass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,qmass,&
735 [ - + ]: 30 : & strten,strtarget,ucvol,vel,vlogv)
736 : :
737 : : !Arguments ------------------------------------
738 : : !scalars
739 : : integer,intent(in) :: nnos,natom
740 : : real(dp),intent(in) :: dtion,ktemp,ucvol,bmass
741 : : real(dp),intent(inout) :: vlogv
742 : : real(dp),intent(out) :: ekin
743 : : type(mttk_type) :: mttk_vars
744 : : !arrays
745 : : real(dp),intent(in) :: amass(natom),strtarget(6),strten(6)
746 : : real(dp),intent(inout) :: vel(3,natom)
747 : : real(dp),intent(in) :: qmass(:)
748 : : integer,intent(in) :: iatfix(:,:)
749 : :
750 : : !Local variables ------------------------------
751 : : !scalars
752 : : integer :: iatom,idir,inos
753 : : real(dp) :: alocal,glogv,gn1kt,gnkt,nfree,odnf,press,prtarget,scale
754 : : logical :: DEBUG=.FALSE.
755 : : !character(len=500) :: message
756 : : !arrays
757 : 30 : real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:)
758 : :
759 : : !***************************************************************************
760 : : !Beginning of executable session
761 : : !***************************************************************************
762 : :
763 [ - + ]: 30 : if(DEBUG) then
764 : 0 : write(std_out,*) ch10,'***',ch10
765 : 0 : write(std_out,*)' isopress : enter '
766 : 0 : write(std_out,*)' strtarget=',strtarget
767 : 0 : write(std_out,*)' bmass=',bmass
768 : 0 : write(std_out,*)' dtion=',dtion
769 : 0 : write(std_out,*)' ktemp=',ktemp
770 : 0 : write(std_out,*)' natom=',natom
771 : 0 : write(std_out,*)' nnos=',nnos
772 : 0 : write(std_out,*)' strten=',strten
773 : 0 : write(std_out,*)' ucvol=',ucvol
774 : : end if
775 : :
776 [ - + + - : 90 : ABI_MALLOC(glogs,(nnos))
+ - + - -
+ - + ]
777 [ - + + - : 60 : ABI_MALLOC(vlogs,(nnos))
- + - + ]
778 [ - + + - : 60 : ABI_MALLOC(xlogs,(nnos))
- + - + ]
779 [ + + ]: 120 : glogs(:)=mttk_vars%glogs(:)
780 [ + + ]: 120 : vlogs(:)=mttk_vars%vlogs(:)
781 [ + + ]: 120 : xlogs(:)=mttk_vars%xlogs(:)
782 : 30 : glogv =mttk_vars%glogv
783 : 30 : scale=one
784 : : !Compute the ionic kinetic energy
785 : 30 : nfree=zero
786 : 30 : ekin=zero
787 [ + + ]: 150 : do iatom=1,natom
788 [ + + ]: 510 : do idir=1,3
789 : : ! Warning : the fixing of atomis is implemented in reduced
790 : : ! coordinates, so that this expression is wrong
791 [ + - ]: 480 : if (iatfix(idir,iatom) == 0) then
792 : 360 : ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2
793 : : ! Counts the degrees of freedom
794 : 360 : nfree=nfree+one
795 : : end if
796 : : end do
797 : : end do
798 : 30 : prtarget=-(strtarget(1)+strtarget(2)+strtarget(3))/three
799 : 30 : press=-(strten(1)+strten(2)+strten(3))/three
800 : 30 : gnkt=nfree*ktemp
801 : 30 : gn1kt=(nfree+one)*ktemp
802 : 30 : odnf=one+three/nfree
803 : : !Update the forces
804 : 30 : glogs(1)=(two*ekin+bmass*vlogv*vlogv-gn1kt)/qmass(1)
805 : 30 : glogv=(odnf*two*ekin+three*(press-prtarget)*ucvol)/bmass
806 : : !Update thermostat velocity
807 : 30 : vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
808 [ + + ]: 90 : do inos=1,nnos-1
809 : 60 : alocal=exp(-dtion/eight*vlogs(nnos+1-inos))
810 : : vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+&
811 : 90 : & dtion/four*glogs(nnos-inos)*alocal
812 : : end do
813 : : !Update dLog(V)/dt
814 : 30 : alocal=exp(-dtion/eight*vlogs(1))
815 : 30 : vlogv=vlogv*alocal**2+dtion/four*glogv*alocal
816 : : !Update the particle velocities
817 : 30 : alocal=exp(-dtion/two*(vlogs(1)+odnf*vlogv))
818 : 30 : scale=scale*alocal
819 : 30 : ekin=ekin*alocal**2
820 : 30 : glogv=(odnf*two*ekin+three*(press-prtarget)*ucvol)/bmass
821 : : !Update the thermostat positions
822 [ + + ]: 120 : do inos=1,nnos
823 : 120 : xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two
824 : : end do
825 : : !Update dLog(V)/dt
826 : 30 : alocal=exp(-dtion/eight*vlogs(1))
827 : 30 : vlogv=vlogv*alocal**2+dtion/four*glogv*alocal
828 : : !Update the forces
829 : 30 : glogs(1)=(two*ekin+bmass*vlogv*vlogv-gn1kt)/qmass(1)
830 : : !Update the thermostat velocities
831 [ + + ]: 90 : do inos=1,nnos-1
832 : 60 : alocal=exp(-dtion/eight*vlogs(inos+1))
833 : 60 : vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal
834 : 90 : glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1)
835 : : end do
836 : 30 : vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
837 [ + + + + ]: 510 : vel(:,:)=vel(:,:)*scale
838 : : !Compute the ionic kinetic energy
839 : 30 : ekin=zero
840 [ + + ]: 150 : do iatom=1,natom
841 [ + + ]: 510 : do idir=1,3
842 : : ! Warning : the fixing of atomis is implemented in reduced
843 : : ! coordinates, so that this expression is wrong
844 [ + - ]: 480 : if (iatfix(idir,iatom) == 0) then
845 : 360 : ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2
846 : : end if
847 : : end do
848 : : end do
849 : : !Compute the thermostat kinetic energy and add it to the ionic one
850 : : !First thermostat
851 : 30 : ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*(nfree+one)*ktemp
852 : : !Other thermostats
853 [ + + ]: 90 : do inos=2,nnos
854 : 90 : ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp
855 : : end do
856 : : !Barostat
857 : 30 : ekin=ekin+half*bmass*vlogv**2+prtarget*ucvol
858 : 30 : ABI_FREE(glogs)
859 : 30 : ABI_FREE(vlogs)
860 : 30 : ABI_FREE(xlogs)
861 : :
862 : : !DEBUG
863 : : !write(std_out,*) 'EKIN',ekin
864 : : !write(std_out,*) 'VLOGV',vlogv
865 : : !write(std_out,*)'ekin added T',half*qmass(:)*vlogs(:)**2,xlogs(:)*(nfree)*ktemp
866 : : !write(std_out,*)'ekin added P',half*bmass*vlogv**2,prtarget*ucvol
867 : : !write(std_out,*)'ekin last',ekin
868 : : !ENDDEBUG
869 : 30 : end subroutine isopress
870 : : !!***
871 : :
872 : : !----------------------------------------------------------------------
873 : :
874 : : !!****f* ABINIT/isostress
875 : : !! NAME
876 : : !! isostress
877 : : !!
878 : : !! FUNCTION
879 : : !! performs one half step on isostress parameters according to Martyna et al.
880 : : !!
881 : : !! INPUTS
882 : : !! amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...)
883 : : !! dtion= ionic time step
884 : : !! isotemp_data
885 : : !! ktemp
886 : : !! press= current pressure of the system
887 : : !! prtarget= target pressure
888 : : !! ucvol= unit cell volume
889 : : !! vel= current velocity
890 : : !!
891 : : !! OUTPUT
892 : : !! Only updates variables
893 : : !!
894 : : !! SIDE EFFECTS
895 : : !! isotemp_data: updates the thermostat parameters (saved variables: bouh !)
896 : : !! vel=update the velocities
897 : : !!
898 : : !! SOURCE
899 : :
900 [ - + ]: 1932 : subroutine isostress(amass,bmass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,&
901 [ - + ]: 1932 : & qmass,strten,strtarget,ucvol,vel)
902 : :
903 : : use m_linalg_interfaces
904 : :
905 : : !Arguments ------------------------------------
906 : : !scalars
907 : : integer,intent(in) :: natom,nnos
908 : : real(dp),intent(in) :: dtion,ktemp,ucvol,bmass
909 : : real(dp),intent(out) :: ekin
910 : : type(mttk_type) :: mttk_vars
911 : : !arrays
912 : : real(dp),intent(in) :: amass(natom),strtarget(6),strten(6),qmass(:)
913 : : real(dp),intent(inout) :: vel(3,natom)
914 : : integer, intent(in) :: iatfix(:,:)
915 : :
916 : : !Local variables ------------------------------
917 : : !scalars
918 : : integer,parameter :: lwork=8
919 : : integer :: iatom,idir,info,inos,jdir
920 : : real(dp) :: akinb,alocal,gn1kt,gnd2kt,nfree,odnf,trvg !,gnkt,scale
921 : : logical :: DEBUG=.FALSE.
922 : : !character(len=500) :: message
923 : : !arrays
924 : : real(dp) :: akin(3,3),expdiag(3),gboxg(3,3),identity(3,3),press(3,3)
925 : : real(dp) :: prtarget(3,3),tvtemp(3,3),uv(3),vboxg(3,3),veig(3),vtemp(3,3)
926 : : real(dp) :: work(lwork)
927 : 1932 : real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:)
928 : :
929 : : !***************************************************************************
930 : : !Beginning of executable session
931 : : !***************************************************************************
932 : :
933 [ - + ]: 1932 : if(DEBUG) then
934 : 0 : write(std_out,*) ch10,'***',ch10
935 : 0 : write(std_out,*)' isostress : enter '
936 : 0 : write(std_out,*)' strtarget=',strtarget
937 : 0 : write(std_out,*)' bmass=',bmass
938 : 0 : write(std_out,*)' dtion=',dtion
939 : 0 : write(std_out,*)' ktemp=',ktemp
940 : 0 : write(std_out,*)' natom=',natom
941 : 0 : write(std_out,*)' nnos=',nnos
942 : 0 : write(std_out,*)' strten=',strten
943 : 0 : write(std_out,*)' ucvol=',ucvol
944 : : end if
945 : :
946 [ - + + - : 5796 : ABI_MALLOC(glogs,(nnos))
+ - + - -
+ - + ]
947 [ - + + - : 3864 : ABI_MALLOC(vlogs,(nnos))
- + - + ]
948 [ - + + - : 3864 : ABI_MALLOC(xlogs,(nnos))
- + - + ]
949 [ + + ]: 3924 : glogs(:)=mttk_vars%glogs(:)
950 [ + + ]: 3924 : vlogs(:)=mttk_vars%vlogs(:)
951 [ + + ]: 3924 : xlogs(:)=mttk_vars%xlogs(:)
952 [ + + + + ]: 25116 : vboxg(:,:)=mttk_vars%vboxg(:,:)
953 : 1932 : identity(:,:)=zero
954 [ + + ]: 7728 : do idir=1,3
955 : 7728 : identity(idir,idir)=one
956 : : end do
957 : :
958 : : !write(std_out,*) 'isostress 02'
959 : : !##########################################################
960 : : !### 03. Compute the ionic kinetic energy
961 : :
962 : 1932 : nfree=zero
963 : 1932 : ekin=zero
964 [ + + ]: 393796 : do iatom=1,natom
965 [ + + ]: 1569388 : do idir=1,3
966 : : ! Warning : the fixing of atomis is implemented in reduced
967 : : ! coordinates, so that this exprtargetion is wrong
968 [ + - ]: 1567456 : if (iatfix(idir,iatom) == 0) then
969 : 1175592 : ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2
970 : : ! Counts the degrees of freedom
971 : 1175592 : nfree=nfree+one
972 : : end if
973 : : end do
974 : : end do
975 : :
976 : 1932 : gn1kt=(nfree+one)*ktemp
977 : 1932 : gnd2kt=(nfree+9)*ktemp
978 : 1932 : odnf=one+three/nfree
979 : 1932 : akin(:,:)=zero
980 [ + + ]: 393796 : do iatom=1,natom
981 [ + + ]: 1569388 : do idir=1,3
982 [ + + ]: 5094232 : do jdir=1,3
983 : : ! Warning : the fixing of atomis is implemented in reduced
984 : : ! coordinates, so that this expression is wrong
985 : 4702368 : akin(idir,jdir)=akin(idir,jdir)+0.5d0*amass(iatom)*vel(idir,iatom)*vel(jdir,iatom)
986 : : end do
987 : : end do
988 : : end do
989 : : akinb=zero
990 [ + + ]: 7728 : do idir=1,3
991 [ + + ]: 25116 : do jdir=1,3
992 : 23184 : akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2
993 : : end do
994 : : end do
995 : : !Compute the pressure: from Voigt to tensor notation+kinetic energy
996 [ + + ]: 7728 : do idir=1,3
997 : 5796 : press(idir,idir)=-strten(idir)
998 : 7728 : prtarget(idir,idir)=-strtarget(idir)
999 : : end do
1000 : 1932 : press(3,2)=-strten(4); press(1,3)=-strten(5); press(2,1)=-strten(6)
1001 : 1932 : prtarget(3,2)=-strtarget(4); prtarget(1,3)=-strtarget(5); prtarget(2,1)=-strtarget(6)
1002 : 1932 : press(2,3)=press(3,2); press(3,1)=press(1,3); press(1,2)=press(2,1)
1003 : 1932 : prtarget(2,3)=prtarget(3,2); prtarget(3,1)=prtarget(1,3); prtarget(1,2)=prtarget(2,1)
1004 : : !Update the forces
1005 : 1932 : glogs(1)=(two*ekin+two*akinb-gnd2kt)/qmass(1)
1006 [ + + + + ]: 25116 : gboxg(:,:)=(two*ekin/nfree*identity(:,:)+two*akin(:,:)+(press(:,:)-prtarget(:,:))*ucvol)/bmass
1007 : : !Update thermostat velocity
1008 [ + - ]: 1932 : if (nnos > 0) vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
1009 [ + + ]: 1992 : do inos=1,nnos-1
1010 : 60 : alocal=exp(-dtion/eight*vlogs(nnos+1-inos))
1011 : : vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+&
1012 : 1992 : & dtion/four*glogs(nnos-inos)*alocal
1013 : : end do
1014 : : !Update box velocity
1015 : 1932 : alocal=exp(-dtion/eight*vlogs(1))
1016 : :
1017 [ - + ]: 1932 : if(DEBUG) then
1018 : 0 : write(std_out,*)' gboxg(:,:)=',gboxg(:,:)
1019 : 0 : write(std_out,*)' vboxg(:,:)=',vboxg(:,:)
1020 : 0 : write(std_out,*)' alocal=',alocal
1021 : : end if
1022 : :
1023 [ + + + + ]: 25116 : vboxg(:,:)=vboxg(:,:)*alocal**2+dtion/four*gboxg(:,:)*alocal
1024 : : !Update the thermostat positions
1025 [ + + ]: 3924 : do inos=1,nnos
1026 : 3924 : xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two
1027 : : end do
1028 : : !Update the particle velocities
1029 : 1932 : trvg=(vboxg(1,1)+vboxg(2,2)+vboxg(3,3))/nfree
1030 [ + + + + ]: 25116 : vtemp(:,:)=vboxg(:,:)+(trvg+vlogs(1))*identity(:,:)
1031 : 1932 : call dsyev('V','U',3,vtemp,3,veig,work,lwork,info)
1032 : : !On exit, we have vtemp=U such that tU vtemp U = veig
1033 [ + + + + ]: 25116 : tvtemp(:,:)=transpose(vtemp)
1034 : :
1035 [ - + ]: 1932 : if(DEBUG) then
1036 : 0 : write(std_out,*)' vboxg(:,:)=',vboxg(:,:)
1037 : 0 : write(std_out,*)' vtemp(:,:)=',vtemp(:,:)
1038 : 0 : write(std_out,*)' veig(:)=',veig(:)
1039 : : end if
1040 : :
1041 : 1932 : expdiag(1)=exp(-veig(1)*dtion/two)
1042 : 1932 : expdiag(2)=exp(-veig(2)*dtion/two)
1043 : 1932 : expdiag(3)=exp(-veig(3)*dtion/two)
1044 : : !if(DEBUG) then
1045 : : ! write(std_out,*)' isostress : expdiag(:)=',expdiag(:) ! Do not remove this line : seems to be needed for g95 compilo
1046 : : !end if
1047 [ + + ]: 393796 : do iatom=1,natom
1048 [ + + + + ]: 5094232 : uv(:)=matmul(tvtemp,vel(:,iatom))
1049 [ + + ]: 1567456 : uv(:)=uv(:)*expdiag(:)
1050 [ + + + + : 6271756 : vel(:,iatom)=matmul(vtemp,uv)
+ + ]
1051 : : end do
1052 : : !Compute the ionic kinetic energy
1053 : 1932 : nfree=zero
1054 : 1932 : ekin=zero
1055 [ + + ]: 393796 : do iatom=1,natom
1056 [ + + ]: 1569388 : do idir=1,3
1057 : : ! Warning : the fixing of atomis is implemented in reduced
1058 : : ! coordinates, so that this expression is wrong
1059 [ + - ]: 1567456 : if (iatfix(idir,iatom) == 0) then
1060 : 1175592 : ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2
1061 [ - + ]: 1175592 : if(DEBUG) then
1062 : 0 : write(std_out,*)'kin',iatom,ekin,vel(idir,iatom)
1063 : : end if
1064 : : ! Counts the degrees of freedom
1065 : 1175592 : nfree=nfree+one
1066 : : end if
1067 : : end do
1068 : : end do
1069 : 1932 : gn1kt=(nfree+one)*ktemp
1070 : 1932 : gnd2kt=(nfree+9)*ktemp
1071 : 1932 : odnf=one+three/nfree
1072 : 1932 : akin(:,:)=zero
1073 [ + + ]: 393796 : do iatom=1,natom
1074 [ + + ]: 1569388 : do idir=1,3
1075 [ + + ]: 5094232 : do jdir=1,3
1076 : : ! Warning : the fixing of atomis is implemented in reduced
1077 : : ! coordinates, so that this expression is wrong
1078 : 4702368 : akin(idir,jdir)=akin(idir,jdir)+0.5d0*amass(iatom)*vel(idir,iatom)*vel(jdir,iatom)
1079 : : end do
1080 : : end do
1081 : : end do
1082 [ + + + + ]: 25116 : gboxg(:,:)=(two*ekin/nfree*identity(:,:)+two*akin(:,:)+(press(:,:)-prtarget(:,:))*ucvol)/bmass
1083 : : !Update box velocity
1084 : : alocal=exp(-dtion/eight*vlogs(1))
1085 [ + + + + ]: 25116 : vboxg(:,:)=vboxg(:,:)*alocal**2+dtion/four*gboxg(:,:)*alocal
1086 : : !Compute the box kinetic energy
1087 : : akinb=zero
1088 [ + + ]: 7728 : do idir=1,3
1089 [ + + ]: 25116 : do jdir=1,3
1090 : 23184 : akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2
1091 : : end do
1092 : : end do
1093 : 1932 : glogs(1)=(two*ekin+two*akinb-gnd2kt)/qmass(1)
1094 : : !Update the thermostat velocities
1095 [ + + ]: 1992 : do inos=1,nnos-1
1096 : 60 : alocal=exp(-dtion/eight*vlogs(inos+1))
1097 : 60 : vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal
1098 : 1992 : glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1)
1099 : : end do
1100 [ + - ]: 1932 : if (nnos > 0) vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
1101 : : !Compute the ionic kinetic energy
1102 : 1932 : ekin=zero
1103 [ + + ]: 393796 : do iatom=1,natom
1104 [ + + ]: 1569388 : do idir=1,3
1105 : : ! Warning : the fixing of atomis is implemented in reduced
1106 : : ! coordinates, so that this expression is wrong
1107 [ + - ]: 1567456 : if (iatfix(idir,iatom) == 0) then
1108 : 1175592 : ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2
1109 : : end if
1110 : : end do
1111 : : end do
1112 : : !Compute the thermostat kinetic energy and add it to the ionic one
1113 : : !First thermostat
1114 : 1932 : ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*(nfree+nine)*ktemp
1115 : : !Other thermostats
1116 [ + + ]: 1992 : do inos=2,nnos
1117 : 1992 : ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp
1118 : : end do
1119 : : !Barostat kinetic energy
1120 : : akinb=zero
1121 [ + + ]: 7728 : do idir=1,3
1122 [ + + ]: 25116 : do jdir=1,3
1123 : 23184 : akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2
1124 : : end do
1125 : : end do
1126 : : !ekin is the invariant minus the potential energy
1127 : 1932 : ekin=ekin+akinb+prtarget(1,1)*ucvol
1128 : :
1129 [ + + + + ]: 25116 : mttk_vars%vboxg(:,:)=vboxg(:,:)
1130 [ + + ]: 3924 : mttk_vars%glogs(:)=glogs(:)
1131 [ + + ]: 3924 : mttk_vars%vlogs(:)=vlogs(:)
1132 [ + + ]: 3924 : mttk_vars%xlogs(:)=xlogs(:)
1133 : 1932 : ABI_FREE(glogs)
1134 : 1932 : ABI_FREE(vlogs)
1135 : 1932 : ABI_FREE(xlogs)
1136 : :
1137 [ - + ]: 1932 : if(DEBUG) then
1138 : : ! write(std_out,*)'ekin added T',half*qmass(:)*vlogs(:)**2,xlogs(:)*(nfree)*ktemp
1139 : : ! write(std_out,*)'ekin added P',akinb,prtarget*ucvol
1140 : : ! write(std_out,*)'ekin last',ekin
1141 : 0 : write(std_out,*) ch10,' exiting from isostress',ch10
1142 : : end if
1143 : :
1144 : 1932 : end subroutine isostress
1145 : : !!***
1146 : :
1147 : : !----------------------------------------------------------------------
1148 : :
1149 : : end module m_pred_isothermal
1150 : : !!***
|
