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