Line data Source code
1 : !!****m* ABINIT/m_pimd
2 : !! NAME
3 : !! m_pimd
4 : !!
5 : !! FUNCTION
6 : !! This module provides several routines and datatypes for the
7 : !! Path-Integral Molecular Dynamics (PIMD) implementation.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2010-2026 ABINIT group (GG,MT)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_pimd
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_dtset
28 : use m_errors
29 : use m_io_tools
30 : use m_random_zbq
31 :
32 : use m_numeric_tools, only : uniformrandom
33 : use m_matrix, only : matr3inv
34 : use m_geometry, only : mkradim
35 :
36 : implicit none
37 :
38 : private
39 :
40 : !public procedures
41 : public :: pimd_init
42 : public :: pimd_nullify
43 : public :: pimd_destroy
44 : public :: pimd_init_qtb
45 : public :: pimd_skip_qtb
46 : public :: pimd_print
47 : public :: pimd_is_restart
48 : public :: pimd_temperature
49 : public :: pimd_initvel
50 : public :: pimd_initvel_from_scratch
51 : public :: pimd_langevin_random
52 : public :: pimd_langevin_random_qtb
53 : public :: pimd_langevin_random_bar
54 : public :: pimd_langevin_random_init
55 : public :: pimd_energies
56 : public :: pimd_forces
57 : public :: pimd_langevin_forces
58 : public :: pimd_nosehoover_forces
59 : public :: pimd_stresses
60 : public :: pimd_diff_stress
61 : public :: pimd_predict_taylor
62 : public :: pimd_predict_verlet
63 : public :: pimd_predict_vel
64 : public :: pimd_predict_velcell
65 : public :: pimd_nosehoover_propagate
66 : public :: pimd_coord_transform
67 : public :: pimd_force_transform
68 : public :: pimd_apply_constraint
69 : public :: pimd_mass_spring
70 : !!***
71 :
72 : !!****t* m_pimd/pimd_type
73 : !! NAME
74 : !! pimd_type
75 : !!
76 : !! FUNCTION
77 : !! Datatype with the variables required to perform PIMD
78 : !!
79 : !! NOTES
80 : !!
81 : !! SOURCE
82 :
83 : type,public :: pimd_type
84 : ! Scalars
85 : integer :: adpimd
86 : integer :: constraint
87 : integer :: irandom
88 : integer :: nnos
89 : integer :: ntypat
90 : integer :: optcell
91 : integer :: pitransform
92 : integer :: traj_unit
93 : integer :: use_qtb
94 : integer :: qtb_file_unit
95 : real(dp) :: adpimd_gamma
96 : real(dp) :: friction
97 : real(dp) :: bmass
98 : real(dp) :: dtion
99 : real(dp) :: frictionbar
100 : ! Arrays
101 : integer ,pointer :: typat(:) ! This pointer is associated with dtset%typat
102 : real(dp),pointer :: amu(:) ! This pointer is associated with dtset%%amu_orig(:,1)
103 : real(dp),pointer :: mdtemp(:) ! This pointer is associated with dtset%mdtemp
104 : real(dp),pointer :: pimass(:) ! This pointer is associated with dtset%pimass
105 : real(dp),pointer :: qmass(:) ! This pointer is associated with dtset%qmass
106 : real(dp),pointer :: strtarget(:) ! This pointer is associated with dtset%strtarget
107 : real(dp),pointer :: wtatcon(:,:,:) ! This pointer is associated with dtset%wtatcon
108 : real(dp),allocatable :: zeta_prev(:,:,:,:)
109 : real(dp),allocatable :: zeta (:,:,:,:)
110 : real(dp),allocatable :: zeta_next(:,:,:,:)
111 : real(dp),allocatable :: dzeta (:,:,:,:)
112 : end type pimd_type
113 : !!***
114 :
115 : CONTAINS !===========================================================
116 : !!***
117 :
118 : !!****f* m_pimd/pimd_init
119 : !! NAME
120 : !! pimd_init
121 : !!
122 : !! FUNCTION
123 : !! Initialize a datastructure of type pimd_type.
124 : !! Open file(s) related to this datastructure.
125 : !!
126 : !! INPUTS
127 : !! dtset <type(dataset_type)>=all input variables in current dataset
128 : !! is_master= TRUE if I am the master process (proc 0)
129 : !!
130 : !! OUTPUT
131 : !!
132 : !! SIDE EFFECTS
133 : !! pimd_param=datastructure of type pimd_type.
134 : !! several parameters for Path-Integral MD.
135 : !!
136 : !! SOURCE
137 :
138 3985 : subroutine pimd_init(dtset,pimd_param,is_master,force_imgmov)
139 :
140 : !Arguments ------------------------------------
141 : !scalars
142 : logical,intent(in) :: is_master
143 : type(dataset_type),target,intent(in) :: dtset
144 : type(pimd_type),intent(inout) :: pimd_param
145 : integer,optional,intent(in) :: force_imgmov
146 : !Local variables-------------------------------
147 : !scalars
148 : integer :: ierr,imgmov
149 : character(len=200) :: msg
150 :
151 : !************************************************************************
152 3985 : if(present(force_imgmov)) then
153 10 : imgmov=force_imgmov
154 : else
155 3975 : imgmov=dtset%imgmov
156 : end if
157 :
158 3985 : call pimd_nullify(pimd_param)
159 :
160 3985 : if((imgmov==9).or.(imgmov==10).or.(imgmov==13))then
161 17 : pimd_param%adpimd = dtset%adpimd
162 17 : pimd_param%constraint = dtset%pimd_constraint
163 17 : pimd_param%irandom = dtset%irandom
164 17 : pimd_param%nnos = dtset%nnos
165 17 : pimd_param%ntypat = dtset%ntypat
166 17 : pimd_param%optcell = dtset%optcell
167 17 : pimd_param%pitransform = dtset%pitransform
168 17 : pimd_param%adpimd_gamma= dtset%adpimd_gamma
169 17 : pimd_param%friction = dtset%friction
170 17 : pimd_param%bmass = dtset%bmass
171 17 : pimd_param%dtion = dtset%dtion
172 17 : pimd_param%frictionbar = dtset%frictionbar
173 17 : pimd_param%mdtemp =>dtset%mdtemp
174 17 : pimd_param%pimass =>dtset%pimass
175 17 : pimd_param%strtarget =>dtset%strtarget
176 17 : pimd_param%amu =>dtset%amu_orig(:,1)
177 17 : pimd_param%qmass =>dtset%qmass
178 17 : pimd_param%typat =>dtset%typat
179 17 : pimd_param%wtatcon =>dtset%wtatcon
180 17 : if(imgmov==10)then
181 0 : pimd_param%use_qtb=1
182 0 : if(is_master)then
183 0 : call pimd_init_qtb(dtset,pimd_param%qtb_file_unit)
184 : end if
185 : end if
186 17 : if(imgmov==13)then
187 15 : ABI_MALLOC(pimd_param%zeta_prev,(3,dtset%natom,dtset%nimage,dtset%nnos))
188 12 : ABI_MALLOC(pimd_param%zeta ,(3,dtset%natom,dtset%nimage,dtset%nnos))
189 12 : ABI_MALLOC(pimd_param%zeta_next,(3,dtset%natom,dtset%nimage,dtset%nnos))
190 12 : ABI_MALLOC(pimd_param%dzeta ,(3,dtset%natom,dtset%nimage,dtset%nnos))
191 828 : pimd_param%zeta_prev=zero
192 828 : pimd_param%zeta =zero
193 828 : pimd_param%zeta_next=zero
194 828 : pimd_param%dzeta =zero
195 : end if
196 17 : if(dtset%useria==37)then
197 0 : ierr=open_file('pimd_traj.dat',msg,newunit=pimd_param%traj_unit,form='unformatted')
198 0 : if (ierr/=0) then
199 0 : ABI_ERROR(msg)
200 : end if
201 : end if
202 : end if
203 :
204 3985 : end subroutine pimd_init
205 : !!***
206 :
207 : !----------------------------------------------------------------------
208 :
209 : !!****f* m_pimd/pimd_nullify
210 : !! NAME
211 : !! pimd_nullify
212 : !!
213 : !! FUNCTION
214 : !! Nullify the content of a datastructure of type pimd_type.
215 : !!
216 : !! INPUTS
217 : !!
218 : !! OUTPUT
219 : !!
220 : !! SIDE EFFECTS
221 : !! pimd_param=datastructure of type pimd_type.
222 : !! several parameters for Path-Integral MD.
223 : !!
224 : !! SOURCE
225 :
226 7960 : subroutine pimd_nullify(pimd_param)
227 :
228 : !Arguments ------------------------------------
229 : !scalars
230 : type(pimd_type),intent(inout) :: pimd_param
231 :
232 : !************************************************************************
233 :
234 7960 : pimd_param%adpimd = 0
235 7960 : pimd_param%constraint = 0
236 7960 : pimd_param%irandom = -1
237 7960 : pimd_param%nnos = -1
238 7960 : pimd_param%ntypat = -1
239 7960 : pimd_param%optcell = -1
240 7960 : pimd_param%pitransform = -1
241 7960 : pimd_param%qtb_file_unit= -1
242 7960 : pimd_param%traj_unit = -1
243 7960 : pimd_param%use_qtb = 0
244 7960 : pimd_param%adpimd_gamma = one
245 7960 : pimd_param%friction = zero
246 7960 : pimd_param%bmass = zero
247 7960 : pimd_param%dtion = zero
248 7960 : pimd_param%frictionbar = zero
249 7960 : nullify(pimd_param%mdtemp)
250 7960 : nullify(pimd_param%pimass)
251 7960 : nullify(pimd_param%strtarget)
252 7960 : nullify(pimd_param%amu)
253 7960 : nullify(pimd_param%qmass)
254 7960 : nullify(pimd_param%typat)
255 7960 : nullify(pimd_param%wtatcon)
256 :
257 7960 : end subroutine pimd_nullify
258 : !!***
259 :
260 : !----------------------------------------------------------------------
261 :
262 : !!****f* m_pimd/pimd_destroy
263 : !! NAME
264 : !! pimd_destroy
265 : !!
266 : !! FUNCTION
267 : !! Destroy the content of a datastructure of type pimd_type.
268 : !! Close open file(s) related to this datastructure.
269 : !!
270 : !! INPUTS
271 : !!
272 : !! OUTPUT
273 : !!
274 : !! SIDE EFFECTS
275 : !! pimd_param=datastructure of type pimd_type.
276 : !! several parameters for PIMD.
277 : !!
278 : !! SOURCE
279 :
280 3975 : subroutine pimd_destroy(pimd_param)
281 :
282 : !Arguments ------------------------------------
283 : !scalars
284 : type(pimd_type),intent(inout) :: pimd_param
285 : !arrays
286 : !Local variables-------------------------------
287 : !scalars
288 : integer :: ierr
289 : character(len=100) :: msg
290 : !arrays
291 :
292 : !************************************************************************
293 :
294 3975 : ABI_SFREE(pimd_param%zeta_prev)
295 3975 : ABI_SFREE(pimd_param%zeta)
296 3975 : ABI_SFREE(pimd_param%zeta_next)
297 3975 : ABI_SFREE(pimd_param%dzeta)
298 :
299 3975 : if (pimd_param%qtb_file_unit>0) then
300 0 : if (is_open(pimd_param%qtb_file_unit)) then
301 0 : ierr=close_unit(pimd_param%qtb_file_unit,msg)
302 : end if
303 : end if
304 :
305 3975 : if (pimd_param%traj_unit>0) then
306 0 : if (is_open(pimd_param%traj_unit)) then
307 0 : ierr=close_unit(pimd_param%traj_unit,msg)
308 : end if
309 : end if
310 :
311 3975 : call pimd_nullify(pimd_param)
312 :
313 3975 : end subroutine pimd_destroy
314 : !!***
315 :
316 : !----------------------------------------------------------------------
317 :
318 : !!****f* m_pimd/pimd_init_qtb
319 : !! NAME
320 : !! pimd_init_qtb
321 : !!
322 : !! FUNCTION
323 : !! Only relevant for PIMD + Quantum Thermal Bath (QTB);
324 : !! Initialize reading of PIQTB random force file.
325 : !! This routine should be called only by master proc.
326 : !!
327 : !! INPUTS
328 : !! dtset <type(dataset_type)>=all input variables in current dataset
329 : !!
330 : !! OUTPUT
331 : !! qtb_file_unit=if a PIQTB_force file exists, return its file unit.
332 : !!
333 : !! SOURCE
334 :
335 0 : subroutine pimd_init_qtb(dtset,qtb_file_unit)
336 :
337 : !Arguments ------------------------------------
338 : !scalars
339 : integer,intent(out) :: qtb_file_unit
340 : type(dataset_type),target,intent(in) :: dtset
341 : !Local variables-------------------------------
342 : !scalars
343 : integer :: ierr,ndof_qtb,ntimimage_qtb,nimage_qtb
344 : real(dp) :: dtion_qtb,mdtemp_qtb
345 : character(len=200) :: msg
346 :
347 : !************************************************************************
348 :
349 : !Try to open PIQTB random force file
350 : ierr=open_file('piqtb_force',msg,newunit=qtb_file_unit,&
351 0 : & form='unformatted',status='old')
352 :
353 : !Read first line of the file
354 0 : read(qtb_file_unit) dtion_qtb,ntimimage_qtb,mdtemp_qtb,nimage_qtb,ndof_qtb
355 :
356 : !Check consistency of the read parameters with ABINIT input file
357 0 : if (abs(dtion_qtb-dtset%dtion)>tol6) then
358 0 : msg='dtion read from piqtb_force file different from dtion in input file!'
359 0 : ABI_ERROR(msg)
360 : end if
361 0 : if (abs(mdtemp_qtb-dtset%mdtemp(2))>tol6) then
362 0 : msg='mdtemp read from piqtb_force file different from mdtemp(2) in input file!'
363 0 : ABI_ERROR(msg)
364 : end if
365 0 : if (ntimimage_qtb<dtset%ntimimage) then
366 0 : msg='ntimimage read from piqtb_force file smaller than ntimimage in input file!'
367 0 : ABI_ERROR(msg)
368 : end if
369 0 : if (nimage_qtb/=dtset%nimage) then
370 0 : msg='nimage read from piqtb_force file different from nimage in input file!'
371 0 : ABI_ERROR(msg)
372 : end if
373 0 : if (ndof_qtb/=3*dtset%natom*dtset%nimage) then
374 0 : msg='Nb of degrees of freedom read from piqtb_force not consistent with input file!'
375 0 : ABI_ERROR(msg)
376 : end if
377 :
378 0 : end subroutine pimd_init_qtb
379 : !!***
380 :
381 : !----------------------------------------------------------------------
382 :
383 : !!****f* m_pimd/pimd_skip_qtb
384 : !! NAME
385 : !! pimd_skip_qtb
386 : !!
387 : !! FUNCTION
388 : !! Only relevant in case of PI-QTB:
389 : !! Skip a line in a QTB random force file
390 : !!
391 : !! INPUTS
392 : !! pimd_param=datastructure of type pimd_type.
393 : !! several parameters for Path-Integral MD.
394 : !!
395 : !! OUTPUT
396 : !!
397 : !! SOURCE
398 :
399 0 : subroutine pimd_skip_qtb(pimd_param)
400 :
401 : !Arguments ------------------------------------
402 : !scalars
403 : type(pimd_type),intent(in) :: pimd_param
404 : !arrays
405 : !Local variables-------------------------------
406 : !scalars
407 : character(len=500) :: msg
408 : !arrays
409 :
410 : !************************************************************************
411 :
412 0 : if (pimd_param%use_qtb==0) return
413 :
414 0 : if (pimd_param%qtb_file_unit<0) then
415 0 : msg='QTB forces file unit should be positive!'
416 0 : ABI_BUG(msg)
417 : end if
418 :
419 : !Skip one line QTB random forces file
420 0 : read(pimd_param%qtb_file_unit)
421 :
422 : end subroutine pimd_skip_qtb
423 : !!***
424 :
425 : !----------------------------------------------------------------------
426 :
427 : !!****f* m_pimd/pimd_is_restart
428 : !! NAME
429 : !! pimd_is_restart
430 : !!
431 : !! FUNCTION
432 : !! Determine whether this is a PIMD restart or not:
433 : !! test on value of velocities and corresponding temperature
434 : !!
435 : !! INPUTS
436 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
437 : !! vel(3,natom,nimage)=velocities for each image of the cell
438 : !!
439 : !! OUTPUT
440 : !! pimd_is_restart=1 if temperature is not zero
441 : !!
442 : !! SIDE EFFECTS
443 : !!
444 : !! SOURCE
445 :
446 8 : function pimd_is_restart(mass,vel,vel_cell)
447 :
448 : !Arguments ------------------------------------
449 : !scalars
450 : integer :: pimd_is_restart
451 : !arrays
452 : real(dp),intent(in) :: mass(:,:),vel(:,:,:)
453 : real(dp),intent(in),optional :: vel_cell(:,:)
454 : !Local variables-------------------------------
455 : !scalars
456 : real(dp),parameter :: zero_temp=tol7
457 : !arrays
458 :
459 : !************************************************************************
460 :
461 8 : pimd_is_restart=0
462 8 : if (pimd_temperature(mass,vel)>zero_temp) pimd_is_restart=1
463 8 : if (present(vel_cell)) then
464 13 : if (maxval(vel_cell)>zero_temp) pimd_is_restart=pimd_is_restart+10
465 : end if
466 :
467 8 : end function pimd_is_restart
468 : !!***
469 :
470 : !----------------------------------------------------------------------
471 :
472 : !!****f* m_pimd/pimd_temperature
473 : !! NAME
474 : !! pimd_temperature
475 : !!
476 : !! FUNCTION
477 : !! Compute temperature from velocities and masses
478 : !!
479 : !! INPUTS
480 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
481 : !! vel(3,natom,nimage)=velocities for each image of the cell
482 : !!
483 : !! OUTPUT
484 : !! pimd_temperature=temperature (from all images of the cell)
485 : !!
486 : !! SIDE EFFECTS
487 : !!
488 : !! SOURCE
489 :
490 355 : function pimd_temperature(mass,vel)
491 :
492 : !Arguments ------------------------------------
493 : !scalars
494 : real(dp) :: pimd_temperature
495 : !arrays
496 : real(dp),intent(in) :: mass(:,:),vel(:,:,:)
497 : !Local variables-------------------------------
498 : !scalars
499 : integer :: iatom,iimage,imass,natom,natom_mass,ndir,nimage,nmass
500 : real(dp) :: v2
501 : character(len=500) :: msg
502 : !arrays
503 :
504 : !************************************************************************
505 :
506 355 : ndir=size(vel,1);natom=size(vel,2);nimage=size(vel,3)
507 355 : natom_mass=size(mass,1);nmass=size(mass,2)
508 355 : if (ndir/=3.or.natom<=0.or.nimage<=0) then
509 0 : msg='Wrong sizes for vel array !'
510 0 : ABI_BUG(msg)
511 : end if
512 355 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=nimage)) then
513 0 : msg='Wrong dimensions for array mass !'
514 0 : ABI_BUG(msg)
515 : end if
516 :
517 355 : v2=zero
518 2019 : do iimage = 1, nimage
519 1664 : imass = min(nmass, iimage)
520 6315 : do iatom = 1, natom
521 18848 : v2 = v2 + sum(vel(:, iatom, iimage)**2) * mass(iatom, imass)
522 : end do
523 : end do
524 355 : pimd_temperature=v2/(dble(3*natom*nimage)*kb_HaK)
525 :
526 355 : end function pimd_temperature
527 : !!***
528 :
529 : !----------------------------------------------------------------------
530 :
531 : !!****f* m_pimd/pimd_print
532 : !! NAME
533 : !! pimd_print
534 : !!
535 : !! FUNCTION
536 : !! Print out results related to PIMD (for given time step)
537 : !!
538 : !! INPUTS
539 : !! constraint=type of constraint eventually applied (on a reaction coordinate)
540 : !! constraint_output(2)=several (real) data to be output when a constraint has been applied
541 : !! eharm=harmonic energy
542 : !! eharm_virial=harmonic energy from virial
543 : !! epot=potential energy
544 : !! forces(3,natom,trotter)=forces on atoms in each cell
545 : !! inertmass(natom)=inertial masses of atoms
546 : !! irestart=1 if this is a restart
547 : !! itimimage=index of time step
548 : !! kt=Thermal energy K_b.T
549 : !! natom=number of atoms
550 : !! optcell=option for cell evolution
551 : !! prtstress=flag for stress tensor printing
552 : !! prtvolimg=printing volume
553 : !! rprim(3,3)=dimensionless real space primitive translations
554 : !! stress(3,3,3)=stress tensor (first dimension corresponds to 3 different estimators)
555 : !! temperature= temperature at present time step
556 : !! traj_unit=flag activating printing of the trajectory in an external file\
557 : !! if >0, indicates the unit number of the trajectory file
558 : !! trotter=Trotter number
559 : !! vel(3,natom,trotter)=velocities of atoms in in each cell
560 : !! vel_cell(3,3)= velocities of cell parameters (time derivative of rprimd)
561 : !! xcart(3,natom,trotter)=cartesian coordinates of atoms in in each cell
562 : !! xred(3,natom,trotter)=reduced coordinates of atoms in in each cell
563 : !! [print_forces]= optional argument; if TRUE print all forces on atoms (default=FALSE)
564 : !!
565 : !! OUTPUT
566 : !! -- only printing --
567 : !!
568 : !! SIDE EFFECTS
569 : !!
570 : !! SOURCE
571 :
572 45 : subroutine pimd_print(constraint,constraint_output,eharm,eharm_virial,epot,&
573 45 : & forces,inertmass,irestart,itimimage,kt,natom,optcell,prtstress,&
574 : & prtvolimg,rprimd,stress,temperature,traj_unit,&
575 45 : & trotter,vel,vel_cell,xcart,xred, &
576 : & print_forces) ! optional argument
577 :
578 : !Arguments ------------------------------------
579 : !scalars
580 : integer,intent(in) :: constraint,irestart,itimimage,natom,optcell
581 : integer,intent(in) :: prtstress,prtvolimg,traj_unit,trotter
582 : logical,optional,intent(in) :: print_forces
583 : real(dp),intent(in) :: eharm,eharm_virial,epot,kt,temperature
584 : !arrays
585 : real(dp),intent(in) :: constraint_output(2)
586 : real(dp),intent(in) :: forces(3,natom,trotter),inertmass(natom)
587 : real(dp),intent(in) :: rprimd(3,3),stress(3,3,3),vel(3,natom,trotter),vel_cell(3,3)
588 : real(dp),intent(in) :: xcart(3,natom,trotter),xred(3,natom,trotter)
589 : !Local variables-------------------------------
590 : !scalars
591 : integer :: iatom,ii,iimage
592 : real(dp) :: mtot
593 : character(len=500) :: msg
594 : !arrays
595 : real(dp) :: acell(3),cdm_cart(3),cdm_red(3),forcetot(3),rprim(3,3)
596 45 : real(dp),allocatable :: centroid(:,:),qudeloc(:)
597 :
598 : !************************************************************************
599 :
600 : !Temperature
601 45 : if(itimimage==1)then
602 8 : msg=ch10
603 8 : if(mod(irestart,10)==0.and.trotter>1) then
604 6 : write(msg,'(2a)') ch10,' This is a PIMD calculation from scratch'
605 2 : else if (mod(irestart,10)==1) then
606 0 : write(msg,'(2a)') ch10,' This is a RESTART calculation'
607 : end if
608 8 : call wrtout(ab_out,msg,'COLL')
609 8 : call wrtout(std_out,msg,'COLL')
610 : ! write(msg,'(a,f12.5,a)') &
611 : !& ' In the initial configuration, the temperature is ',temperature1,' K'
612 : ! call wrtout(ab_out,msg,'COLL')
613 : ! call wrtout(std_out,msg,'COLL')
614 : end if
615 45 : if (trotter>1) then
616 35 : write(msg,'(2a,i5,a,f12.5,a)') ch10,&
617 70 : & ' At PIMD time step ',itimimage,', the temperature is',temperature,' K'
618 : else
619 10 : write(msg,'(2a,i5,a,f12.5,a)') ch10,&
620 20 : & ' At time step ',itimimage,', the temperature is',temperature,' K'
621 : end if
622 45 : call wrtout(ab_out,msg,'COLL')
623 45 : call wrtout(std_out,msg,'COLL')
624 :
625 : !Energies
626 45 : if (trotter>1) then
627 35 : write(msg,'(4a,f18.9,a,a,a,f18.9,a,a)') ch10,&
628 35 : & ' Energy:',ch10, &
629 35 : & ' Internal energy (PRIMITIVE estimator) =',onehalf*dble(natom*trotter)*kt-eharm+epot ,' Ha',ch10, &
630 70 : & ' Internal energy (VIRIAL estimator) =',onehalf*dble(natom)*kt+eharm_virial+epot,' Ha',ch10
631 : else
632 10 : write(msg,'(2a,f18.9,a)') ch10,&
633 20 : & ' Internal energy =',onehalf*dble(natom)*temperature*kb_HaK+epot ,' Ha'
634 : end if
635 45 : call wrtout(ab_out,msg,'COLL')
636 45 : call wrtout(std_out,msg,'COLL')
637 :
638 : !Forces
639 45 : if (present(print_forces)) then
640 0 : if (print_forces) then
641 0 : write(msg,'(2a)') ch10,' Forces on atoms:'
642 0 : call wrtout(std_out,msg,'COLL')
643 0 : call wrtout(ab_out,msg,'COLL')
644 0 : do iimage=1,trotter
645 0 : select case(iimage)
646 : case(1)
647 0 : write(msg,'(a)') ' fcart'
648 : case(2,3,4,5,6,7,8,9)
649 0 : write(msg,'(a,i1,a)') ' fcart_',iimage,'img'
650 : case(10:99)
651 0 : write(msg,'(a,i2,a)') ' fcart_',iimage,'img'
652 : case default
653 0 : write(msg,'(a,i3,a)') ' fcart_',iimage,'img'
654 : end select
655 0 : call wrtout(std_out,msg,'COLL')
656 0 : call wrtout(ab_out,msg,'COLL')
657 0 : do iatom=1,natom
658 0 : write(msg,'(3f18.10)') forces(1:3,iatom,iimage)
659 0 : call wrtout(std_out,msg,'COLL')
660 0 : call wrtout(ab_out,msg,'COLL')
661 : end do
662 : end do
663 : end if
664 : end if
665 :
666 : !Stress tensor and pressure
667 45 : write(msg,'(2a,3(2a,3f18.9))') ch10,&
668 45 : & ' Stress tensor from PRIMITIVE estimator (Ha/Bohr^3):',ch10, &
669 45 : & ' ',stress(1,1,1),stress(1,1,2),stress(1,1,3),ch10, &
670 45 : & ' ',stress(1,2,1),stress(1,2,2),stress(1,2,3),ch10, &
671 90 : & ' ',stress(1,3,1),stress(1,3,2),stress(1,3,3)
672 45 : if (prtstress==1) then
673 40 : call wrtout(ab_out,msg,'COLL')
674 : end if
675 45 : call wrtout(std_out,msg,'COLL')
676 45 : write(msg,'(a,f18.9,a)') ' Pressure (primitive estimator) =', &
677 90 : & -third*(stress(1,1,1)+stress(1,2,2)+stress(1,3,3))*HaBohr3_GPa,' GPa'
678 45 : call wrtout(ab_out,msg,'COLL')
679 45 : call wrtout(std_out,msg,'COLL')
680 :
681 : !Data related to constraint eventually applied
682 45 : if (constraint/=0) then
683 5 : if (constraint==1) write(msg,'(2a)') ch10,' Blue Moon Ensemble method is activated:'
684 5 : call wrtout(ab_out,msg,'COLL')
685 5 : call wrtout(std_out,msg,'COLL')
686 : write(msg,'(a,f18.10,2a,f18.10)') &
687 5 : & ' - Reaction coordinate =',constraint_output(1),ch10,&
688 10 : & ' - Instantaneous force on the reaction coord. =',constraint_output(2)
689 5 : call wrtout(ab_out,msg,'COLL')
690 5 : call wrtout(std_out,msg,'COLL')
691 : end if
692 :
693 : !Total force
694 45 : if (prtvolimg<=1) then
695 40 : forcetot=zero
696 230 : do iimage=1,trotter
697 650 : do iatom=1,natom
698 1870 : do ii=1,3
699 1680 : forcetot(ii)=forcetot(ii)+forces(ii,iatom,iimage)
700 : end do
701 : end do
702 : end do
703 40 : write(msg,'(2a,3f18.10)') ch10,' Total force=',forcetot(1:3)
704 40 : call wrtout(std_out,msg,'COLL')
705 : end if
706 :
707 : !position of mass center
708 45 : if (prtvolimg<=1) then
709 40 : mtot=zero;cdm_cart=zero;cdm_red=zero
710 230 : do iimage=1,trotter
711 650 : do iatom=1,natom
712 1680 : cdm_cart(:)=cdm_cart(:)+inertmass(iatom)*xcart(:,iatom,iimage)
713 1680 : cdm_red (:)=cdm_red (:)+inertmass(iatom)*xred (:,iatom,iimage)
714 610 : mtot=mtot+inertmass(iatom)
715 : end do
716 : end do
717 160 : cdm_cart=cdm_cart/mtot
718 160 : cdm_red =cdm_red /mtot
719 40 : write(msg,'(3a,3x,3f18.10,3a,3x,3f18.10)') ch10,&
720 40 : & ' Center of mass, in cartes. coordinates :',ch10,cdm_cart(:),ch10,&
721 80 : & ' Center of mass, in reduced coordinates :',ch10,cdm_red(:)
722 40 : call wrtout(std_out,msg,'COLL')
723 40 : call wrtout(ab_out,msg,'COLL')
724 : end if
725 :
726 : !Positions
727 45 : write(msg,'(2a)') ch10,' Atomic positions:'
728 45 : call wrtout(std_out,msg,'COLL')
729 45 : call wrtout(ab_out,msg,'COLL')
730 260 : do iimage=1,trotter
731 45 : select case(iimage)
732 : case(1)
733 45 : write(msg,'(a)') ' xred'
734 : case(2,3,4,5,6,7,8,9)
735 170 : write(msg,'(a,i1,a)') ' xred_',iimage,'img'
736 : case(10:99)
737 0 : write(msg,'(a,i2,a)') ' xred_',iimage,'img'
738 : case default
739 215 : write(msg,'(a,i3,a)') ' xred_',iimage,'img'
740 : end select
741 215 : call wrtout(std_out,msg,'COLL')
742 215 : call wrtout(ab_out,msg,'COLL')
743 215 : if (traj_unit>0) then
744 0 : call wrtout(traj_unit,msg,'COLL')
745 : end if
746 830 : do iatom=1,natom
747 570 : write(msg,'(3f18.10)') xred(1:3,iatom,iimage)
748 570 : call wrtout(std_out,msg,'COLL')
749 570 : call wrtout(ab_out,msg,'COLL')
750 785 : if (traj_unit>0) then
751 0 : call wrtout(traj_unit,msg,'COLL')
752 : end if
753 : end do
754 : end do
755 :
756 : !Velocities
757 45 : write(msg,'(2a)') ch10,' Velocities:'
758 45 : call wrtout(std_out,msg,'COLL')
759 45 : call wrtout(ab_out,msg,'COLL')
760 260 : do iimage=1,trotter
761 45 : select case(iimage)
762 : case(1)
763 45 : write(msg,'(a)') ' vel'
764 : case(2,3,4,5,6,7,8,9)
765 170 : write(msg,'(a,i1,a)') ' vel_',iimage,'img'
766 : case(10:99)
767 0 : write(msg,'(a,i2,a)') ' vel_',iimage,'img'
768 : case default
769 215 : write(msg,'(a,i3,a)') ' vel_',iimage,'img'
770 : end select
771 215 : call wrtout(std_out,msg,'COLL')
772 215 : call wrtout(ab_out,msg,'COLL')
773 215 : if (traj_unit>0) then
774 0 : call wrtout(traj_unit,msg,'COLL')
775 : end if
776 830 : do iatom=1,natom
777 570 : write(msg,'(3f18.10)') vel(1:3,iatom,iimage)
778 570 : call wrtout(std_out,msg,'COLL')
779 570 : call wrtout(ab_out,msg,'COLL')
780 785 : if (traj_unit>0) then
781 0 : call wrtout(traj_unit,msg,'COLL')
782 : end if
783 : end do
784 : end do
785 :
786 45 : if (optcell>0) then
787 :
788 5 : call mkradim(acell,rprim,rprimd)
789 :
790 : ! Time derivative of rprimd
791 5 : write(msg,'(2a)') ch10,' Time derivative of rprimd:'
792 5 : call wrtout(std_out,msg,'COLL')
793 5 : call wrtout(ab_out,msg,'COLL')
794 5 : write(msg,'(2a,3(3f18.10))') ' vel_cell',ch10,vel_cell(:,:)
795 5 : call wrtout(std_out,msg,'COLL')
796 5 : call wrtout(ab_out,msg,'COLL')
797 5 : if (traj_unit>0) then
798 0 : call wrtout(traj_unit,msg,'COLL')
799 : end if
800 :
801 : ! rprimd
802 5 : write(msg,'(2a)') ch10,' Cell parameters:'
803 5 : call wrtout(std_out,msg,'COLL')
804 5 : call wrtout(ab_out,msg,'COLL')
805 5 : write(msg,'(2a,3(3f18.10),3a,3f18.10)') ' rprim',ch10,rprim(:,:),ch10,&
806 10 : & ' acell',ch10,acell(:)
807 5 : call wrtout(std_out,msg,'COLL')
808 5 : call wrtout(ab_out,msg,'COLL')
809 5 : if (traj_unit>0) then
810 0 : call wrtout(traj_unit,msg,'COLL')
811 : end if
812 :
813 : end if
814 :
815 : !Centroids and wave-packet spatial spreads
816 135 : ABI_MALLOC(centroid,(3,natom))
817 135 : ABI_MALLOC(qudeloc,(natom))
818 795 : centroid=zero;qudeloc=zero
819 260 : do iimage=1,trotter
820 830 : do iatom=1,natom
821 2495 : do ii=1,3
822 2280 : centroid(ii,iatom)=centroid(ii,iatom)+xcart(ii,iatom,iimage)
823 : end do
824 : end do
825 : end do
826 645 : centroid=centroid/dble(trotter)
827 260 : do iimage=1,trotter
828 830 : do iatom=1,natom
829 2495 : do ii=1,3
830 2280 : qudeloc(iatom)=qudeloc(iatom)+((xcart(ii,iatom,iimage)-centroid(ii,iatom))**2)
831 : end do
832 : end do
833 : end do
834 195 : qudeloc(:)=sqrt(qudeloc(:)/dble(trotter))
835 45 : write(msg,'(4a)') ch10,' Centroids and wave-packet spatial spreads (cart. coord.):',ch10,&
836 90 : & ' iat centroid_x centroid_y centroid_z spatial_spread'
837 45 : call wrtout(std_out,msg,'COLL')
838 195 : do iatom=1,natom
839 150 : write(msg,'(i4,4f18.10)') iatom,centroid(1:3,iatom),qudeloc(iatom)
840 195 : call wrtout(std_out,msg,'COLL')
841 : end do
842 45 : ABI_FREE(centroid)
843 45 : ABI_FREE(qudeloc)
844 :
845 : !Fake statement
846 : return;ii=prtvolimg
847 :
848 : end subroutine pimd_print
849 : !!***
850 :
851 : !----------------------------------------------------------------------
852 :
853 : !!****f* m_pimd/pimd_initvel
854 : !! NAME
855 : !! pimd_initvel
856 : !!
857 : !! FUNCTION
858 : !! Initialize velocities for PIMD with a gaussian distribution
859 : !! fixing the center of mass
860 : !! and eventually applying a constraint on atomic positions
861 : !!
862 : !! INPUTS
863 : !! constraint=type of constraint to be applied
864 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
865 : !! natom=number of atoms
866 : !! temperature=temperature used to define velocities
867 : !! trotter=Trotter number
868 : !! wtatcon(3,natom)=weights for atomic constraints
869 : !!
870 : !! OUTPUT
871 : !! vel(3,natom,trotter)=velocities of atoms in in each cell
872 : !!
873 : !! SIDE EFFECTS
874 : !! iseed=seed for random number generator
875 : !!
876 : !! SOURCE
877 :
878 8 : subroutine pimd_initvel(iseed,mass,natom,temperature,trotter,vel,constraint,wtatcon)
879 :
880 : !Arguments ------------------------------------
881 : !scalars
882 : integer,intent(in) :: constraint,natom,trotter
883 : integer,intent(inout) :: iseed
884 : real(dp),intent(in) :: temperature
885 : !arrays
886 : real(dp),intent(in) :: mass(:,:),wtatcon(3,natom)
887 : real(dp),intent(out) :: vel(3,natom,trotter)
888 : !Local variables-------------------------------
889 : !scalars
890 : integer :: iatom,ii,iimage,imass,natom_mass,nmass
891 : real(dp) :: mtot,rescale_vel
892 : character(len=500) :: msg
893 : !arrays
894 : real(dp) :: mvini(3)
895 :
896 : !************************************************************************
897 :
898 8 : natom_mass=size(mass,1);nmass=size(mass,2)
899 8 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
900 0 : msg='Wrong dimensions for array mass !'
901 0 : ABI_BUG(msg)
902 : end if
903 :
904 : !Compute total mass (of non-constrained atoms)
905 8 : if (constraint==0) then
906 100 : mtot=sum(mass(1:natom,1:nmass))
907 : else
908 : mtot=zero
909 0 : do iatom=1,natom
910 0 : if (all(abs(wtatcon(:,iatom))<tol8)) mtot=mtot+sum(mass(iatom,1:nmass))
911 : end do
912 : end if
913 8 : if (nmass==1) mtot=mtot*dble(trotter)
914 :
915 : !Initialize randomly the velocities
916 46 : do iimage=1,trotter
917 38 : imass=min(nmass,iimage)
918 130 : do iatom=1,natom
919 374 : do ii=1,3
920 252 : vel(ii,iatom,iimage)=sqrt(kb_HaK*temperature/mass(iatom,imass))*cos(two_pi*uniformrandom(iseed))
921 336 : vel(ii,iatom,iimage)=vel(ii,iatom,iimage)*sqrt(-two*log(uniformrandom(iseed)))
922 : end do
923 : end do
924 : end do
925 :
926 : !Cancel velocities of constrained atoms
927 8 : if (constraint/=0) then
928 0 : do iimage=1,trotter
929 0 : do iatom=1,natom
930 0 : if (any(abs(wtatcon(:,iatom))>=tol8)) vel(:,iatom,iimage)=zero
931 : end do
932 : end do
933 : end if
934 :
935 : !Make sure that the (sum of m_i v_i) at step zero is zero
936 8 : mvini=zero
937 46 : do iimage=1,trotter
938 38 : imass=min(nmass,iimage)
939 130 : do iatom=1,natom
940 374 : do ii=1,3
941 336 : mvini(ii)=mvini(ii)+mass(iatom,imass)*vel(ii,iatom,iimage)
942 : end do
943 : end do
944 : end do
945 8 : if (constraint==0) then
946 46 : do iimage=1,trotter
947 130 : do iatom=1,natom
948 374 : do ii=1,3
949 336 : vel(ii,iatom,iimage)=vel(ii,iatom,iimage)-(mvini(ii)/mtot)
950 : end do
951 : end do
952 : end do
953 : else
954 0 : do iimage=1,trotter
955 0 : do iatom=1,natom
956 0 : if (all(abs(wtatcon(:,iatom))<tol8)) vel(:,iatom,iimage)=vel(:,iatom,iimage)-(mvini(:)/mtot)
957 : end do
958 : end do
959 : end if
960 :
961 : !Now rescale the velocities to give the exact temperature
962 8 : rescale_vel=sqrt(temperature/pimd_temperature(mass,vel))
963 382 : vel(:,:,:)=vel(:,:,:)*rescale_vel
964 :
965 8 : end subroutine pimd_initvel
966 : !!***
967 :
968 : !----------------------------------------------------------------------
969 :
970 : !!****f* m_pimd/pimd_initvel_from_scratch
971 : !! NAME
972 : !! pimd_initvel_from_scratch
973 : !!
974 : !! FUNCTION
975 : !! Initialize velocities for PIMD with a gaussian distribution
976 : !! at start of the calculation (with initial temperature)
977 : !!
978 : !! INPUTS
979 : !! pimd_temperature=temperature (from all images of the cell)
980 : !! natom=number of atoms
981 : !! trotter=Trotter number
982 : !! vel_cell(3,3)=cell velocities
983 : !!
984 : !! OUTPUT
985 : !! vel(3,natom,nimage)=velocities for each image of the cell
986 : !!
987 : !! SOURCE
988 :
989 0 : subroutine pimd_initvel_from_scratch(pimd_param,natom,trotter,vel,vel_cell)
990 :
991 : !Arguments ------------------------------------
992 : !scalars
993 : integer,intent(in) :: natom,trotter
994 : real(dp),intent(in) :: vel_cell(3,3)
995 : real(dp),intent(out) :: vel(3,natom,trotter)
996 : type(pimd_type),intent(inout) :: pimd_param
997 : !Local variables-------------------------------
998 : !scalars
999 : integer :: iseed=-5,iimage,irestart,ndof,zeroforce
1000 : real(dp) :: initemp,rescale_temp
1001 : !arrays
1002 0 : real(dp),allocatable :: inertmass(:,:)
1003 :
1004 : !************************************************************************
1005 :
1006 0 : ABI_MALLOC(inertmass,(natom,trotter))
1007 0 : do iimage=1,trotter
1008 0 : inertmass(1:natom,iimage)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
1009 : end do
1010 :
1011 0 : irestart=pimd_is_restart(inertmass,vel,vel_cell)
1012 :
1013 : !Initialize derivatives
1014 0 : if (mod(irestart,10)==0) then
1015 0 : ndof=3*natom*trotter
1016 0 : rescale_temp=one ; if(zeroforce==1) rescale_temp=dble(ndof)/dble(ndof-3)
1017 : zeroforce=1 ; if(pimd_param%pitransform==1.or.pimd_param%pitransform==2) zeroforce=0
1018 0 : if(pimd_param%constraint==1) zeroforce=0
1019 0 : initemp=pimd_param%mdtemp(1)/rescale_temp
1020 0 : call pimd_initvel(iseed,inertmass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
1021 : end if
1022 :
1023 0 : ABI_FREE(inertmass)
1024 :
1025 0 : end subroutine pimd_initvel_from_scratch
1026 : !!***
1027 :
1028 : !----------------------------------------------------------------------
1029 :
1030 : !!****f* m_pimd/pimd_langevin_random
1031 : !! NAME
1032 : !! pimd_langevin_random
1033 : !!
1034 : !! FUNCTION
1035 : !! Generate a set of random numbers to be used for PIMD Langevin algorithm
1036 : !!
1037 : !! INPUTS
1038 : !! irandom=option for random number generator:
1039 : !! 1:uniform random routine provided within Abinit package
1040 : !! 2:Fortran 90 random number generator
1041 : !! 3:ZBQLU01 non deterministic random number generator
1042 : !! langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter)
1043 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
1044 : !! natom=number of atoms
1045 : !! trotter=Trotter number
1046 : !! zeroforce=flag; if 1 keep sum of forces equal to zero
1047 : !!
1048 : !! OUTPUT
1049 : !! alea(3,natom,trotter)=set of random numbers
1050 : !!
1051 : !! SIDE EFFECTS
1052 : !! iseed=seed for random number generator (used only if irandom=1)
1053 : !!
1054 : !! SOURCE
1055 :
1056 30 : subroutine pimd_langevin_random(alea,irandom,iseed,langev,mass,natom,trotter,zeroforce)
1057 :
1058 : !Arguments ------------------------------------
1059 : !scalars
1060 : integer,intent(in) :: irandom,natom,trotter,zeroforce
1061 : integer,intent(inout) :: iseed
1062 : !arrays
1063 : real(dp),intent(in) :: langev(:,:),mass(:,:)
1064 : real(dp),intent(out) :: alea(3,natom,trotter)
1065 : !Local variables-------------------------------
1066 : !scalars
1067 : integer :: iatom,ii,iimage,imass,nmass,natom_mass
1068 : real(dp) :: mtot,r1,r2
1069 : character(len=500) :: msg
1070 : !arrays
1071 : real(dp) :: total(3)
1072 :
1073 : !************************************************************************
1074 :
1075 30 : natom_mass=size(mass,1);nmass=size(mass,2)
1076 30 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1077 0 : msg='Wrong dimensions for array mass !'
1078 0 : ABI_BUG(msg)
1079 : end if
1080 :
1081 : mtot=sum(mass(1:natom,1:nmass))
1082 30 : if (nmass==1) mtot=mtot*trotter
1083 :
1084 : !Draw random numbers
1085 155 : do iimage=1,trotter
1086 545 : do iatom=1,natom
1087 1685 : do ii=1,3
1088 2340 : select case(irandom)
1089 : case(1)
1090 1170 : r1=uniformrandom(iseed)
1091 1170 : r2=uniformrandom(iseed)
1092 : case(2)
1093 0 : call random_number(r1)
1094 0 : call random_number(r2)
1095 : case(3)
1096 0 : r1=ZBQLU01(zero)
1097 1170 : r2=ZBQLU01(zero)
1098 : end select
1099 1560 : alea(ii,iatom,iimage)= cos(two_pi*r1)*sqrt(-log(r2)*two)
1100 : end do
1101 : end do
1102 : end do
1103 :
1104 : !Make sure that the sum of random forces is zero
1105 30 : if(zeroforce==1)then
1106 15 : total=zero
1107 15 : mtot=zero
1108 55 : do iimage=1,trotter
1109 40 : imass=min(nmass,iimage)
1110 175 : do iatom=1,natom
1111 160 : mtot=mtot+mass(iatom,imass)
1112 : end do
1113 : end do
1114 55 : do iimage=1,trotter
1115 40 : imass=min(nmass,iimage)
1116 175 : do iatom=1,natom
1117 520 : do ii=1,3
1118 480 : total(ii)=total(ii)+langev(iatom,imass)*alea(ii,iatom,iimage)
1119 : end do
1120 : end do
1121 : end do
1122 55 : do iimage=1,trotter
1123 40 : imass=min(nmass,iimage)
1124 175 : do iatom=1,natom
1125 520 : do ii=1,3
1126 : alea(ii,iatom,iimage)= alea(ii,iatom,iimage)- &
1127 480 : & (total(ii)*mass(iatom,imass))/(langev(iatom,imass)*mtot)
1128 : end do
1129 : end do
1130 : end do
1131 : ! now random forces have been rescaled so that their sum is zero
1132 : end if
1133 :
1134 30 : end subroutine pimd_langevin_random
1135 : !!***
1136 :
1137 : !----------------------------------------------------------------------
1138 :
1139 : !!****f* m_pimd/pimd_langevin_random_qtb
1140 : !! NAME
1141 : !! pimd_langevin_random_qtb
1142 : !!
1143 : !! FUNCTION
1144 : !! Read a set of random forces (atm units) to be used for PIMD QTB algorithm
1145 : !!
1146 : !! INPUTS
1147 : !! langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter)
1148 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
1149 : !! natom=number of atoms
1150 : !! qtb_file_unit= ramdom forces file unit
1151 : !! trotter=Trotter number
1152 : !! zeroforce=flag; if 1 keep sum of forces equal to zero
1153 : !!
1154 : !! OUTPUT
1155 : !! alea(3,natom,trotter)=set of random forces
1156 : !!
1157 : !! SOURCE
1158 :
1159 0 : subroutine pimd_langevin_random_qtb(alea,langev,mass,natom,qtb_file_unit,trotter,zeroforce)
1160 :
1161 : !Arguments ------------------------------------
1162 : !scalars
1163 : integer,intent(in) :: natom,qtb_file_unit,trotter,zeroforce
1164 : !arrays
1165 : real(dp),intent(in) :: langev(:,:),mass(:,:)
1166 : real(dp),intent(out) :: alea(3,natom,trotter)
1167 : !Local variables-------------------------------
1168 : !scalars
1169 : integer :: iatom,ii,iimage,imass,nmass,natom_mass
1170 : real(dp) :: mtot
1171 : character(len=500) :: msg
1172 : !arrays
1173 0 : real(sp) :: alea_sp(3,natom,trotter)
1174 : real(dp) :: total(3)
1175 :
1176 : !************************************************************************
1177 :
1178 0 : if (qtb_file_unit<0) then
1179 0 : msg='QTB forces file unit should be positive!'
1180 0 : ABI_BUG(msg)
1181 : end if
1182 0 : natom_mass=size(mass,1);nmass=size(mass,2)
1183 0 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1184 0 : msg='Wrong dimensions for array mass !'
1185 0 : ABI_BUG(msg)
1186 : end if
1187 :
1188 : !Read QTB random forces
1189 0 : read(qtb_file_unit) alea_sp(1:3,1:natom,1:trotter)
1190 0 : alea(:,:,:)=dble(alea_sp(:,:,:))
1191 :
1192 : !Make sure that the sum of random forces is zero
1193 0 : if(zeroforce==1)then
1194 0 : total=zero
1195 0 : mtot=zero
1196 0 : do iimage=1,trotter
1197 0 : imass=min(nmass,iimage)
1198 0 : do iatom=1,natom
1199 0 : mtot=mtot+mass(iatom,imass)
1200 : end do
1201 : end do
1202 0 : do iimage=1,trotter
1203 0 : imass=min(nmass,iimage)
1204 0 : do iatom=1,natom
1205 0 : do ii=1,3
1206 0 : total(ii)=total(ii)+langev(iatom,imass)*alea(ii,iatom,iimage)
1207 : end do
1208 : end do
1209 : end do
1210 0 : do iimage=1,trotter
1211 0 : imass=min(nmass,iimage)
1212 0 : do iatom=1,natom
1213 0 : do ii=1,3
1214 : alea(ii,iatom,iimage)= alea(ii,iatom,iimage)- &
1215 0 : & (total(ii)*mass(iatom,imass))/(langev(iatom,imass)*mtot)
1216 : end do
1217 : end do
1218 : end do
1219 : ! now random forces have been rescaled so that their sum is zero
1220 : end if
1221 :
1222 0 : end subroutine pimd_langevin_random_qtb
1223 : !!***
1224 :
1225 : !----------------------------------------------------------------------
1226 :
1227 : !!****f* m_pimd/pimd_langevin_random_bar
1228 : !! NAME
1229 : !! pimd_langevin_random_bar
1230 : !!
1231 : !! FUNCTION
1232 : !! Generate a set of random numbers to be used for the barostat of PIMD Langevin algorithm
1233 : !!
1234 : !! INPUTS
1235 : !! irandom=option for random number generator:
1236 : !! 1:uniform random routine provided within Abinit package
1237 : !! 2:Fortran 90 random number generator
1238 : !! 3:ZBQLU01 non deterministic random number generator
1239 : !!
1240 : !! OUTPUT
1241 : !! alea_bar(3,3)=set of random numbers
1242 : !!
1243 : !! SIDE EFFECTS
1244 : !! iseed=seed for random number generator (used only if irandom=1)
1245 : !!
1246 : !! SOURCE
1247 :
1248 5 : subroutine pimd_langevin_random_bar(alea_bar,irandom,iseed)
1249 :
1250 : !Arguments ------------------------------------
1251 : !scalars
1252 : integer,intent(in) :: irandom
1253 : integer,intent(inout) :: iseed
1254 : !arrays
1255 : real(dp),intent(out) :: alea_bar(3,3)
1256 : !Local variables-------------------------------
1257 : !scalars
1258 : integer :: ii,jj
1259 : real(dp) :: r1,r2
1260 : !arrays
1261 :
1262 : !************************************************************************
1263 :
1264 : !Draw random numbers
1265 20 : do ii=1,3
1266 65 : do jj=1,3
1267 90 : select case(irandom)
1268 : case(1)
1269 45 : r1=uniformrandom(iseed)
1270 45 : r2=uniformrandom(iseed)
1271 : case(2)
1272 0 : call random_number(r1)
1273 0 : call random_number(r2)
1274 : case(3)
1275 0 : r1=ZBQLU01(zero)
1276 45 : r2=ZBQLU01(zero)
1277 : end select
1278 60 : alea_bar(ii,jj)= cos(two_pi*r1)*sqrt(-log(r2)*two)
1279 : end do
1280 : end do
1281 :
1282 : !Symmetrize
1283 5 : alea_bar(1,2)=half*(alea_bar(1,2)+alea_bar(2,1))
1284 5 : alea_bar(1,3)=half*(alea_bar(1,3)+alea_bar(3,1))
1285 5 : alea_bar(2,3)=half*(alea_bar(2,3)+alea_bar(3,2))
1286 5 : alea_bar(2,1)=alea_bar(1,2)
1287 5 : alea_bar(3,1)=alea_bar(1,3)
1288 5 : alea_bar(3,2)=alea_bar(2,3)
1289 :
1290 5 : end subroutine pimd_langevin_random_bar
1291 : !!***
1292 :
1293 : !----------------------------------------------------------------------
1294 :
1295 : !!****f* m_pimd/pimd_langevin_random_init
1296 : !! NAME
1297 : !! pimd_langevin_random_init
1298 : !!
1299 : !! FUNCTION
1300 : !! Initialize random number generator to be used for PIMD Langevin algorithm
1301 : !!
1302 : !! INPUTS
1303 : !! irandom=option for random number generator:
1304 : !! 1:uniform random routine provided within Abinit package
1305 : !! 2:Fortran 90 random number generator
1306 : !! 3:ZBQLU01 non deterministic random number generator
1307 : !!
1308 : !! OUTPUT
1309 : !!
1310 : !! SIDE EFFECTS
1311 : !! iseed=seed for random number generator (used only if irandom=1)
1312 : !!
1313 : !! SOURCE
1314 :
1315 5 : subroutine pimd_langevin_random_init(irandom,iseed)
1316 :
1317 : !Arguments ------------------------------------
1318 : !scalars
1319 : integer,intent(in) :: irandom
1320 : integer,intent(inout) :: iseed
1321 :
1322 : !************************************************************************
1323 :
1324 5 : if (irandom==3) then
1325 0 : call ZBQLINI(0)
1326 : end if
1327 :
1328 : !Fake statement
1329 : return;if (.false.) iseed=zero
1330 :
1331 : end subroutine pimd_langevin_random_init
1332 : !!***
1333 :
1334 : !----------------------------------------------------------------------
1335 :
1336 : !!****f* m_pimd/pimd_energies
1337 : !! NAME
1338 : !! pimd_energies
1339 : !!
1340 : !! FUNCTION
1341 : !! In the case od PIMD, compute the several contribution to total energy
1342 : !!
1343 : !! INPUTS
1344 : !! etotal_img(trotter)= energy (from DFT) for each cell
1345 : !! forces(3,natom,trotter)=forces (from DFT) on atoms in each cell
1346 : !! natom=number of atoms
1347 : !! spring(natom)=spring constants in the primitive scheme
1348 : !! trotter=Trotter number
1349 : !! xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
1350 : !!
1351 : !! OUTPUT
1352 : !! eharm =harmonic energy
1353 : !! eharm_virial=harmonic energy from virial estimator
1354 : !! epot =potential energy
1355 : !!
1356 : !! SIDE EFFECTS
1357 : !!
1358 : !! SOURCE
1359 :
1360 45 : subroutine pimd_energies(eharm,eharm_virial,epot,etotal_img,forces,natom,spring,trotter,xcart)
1361 :
1362 : !Arguments ------------------------------------
1363 : !scalars
1364 : integer,intent(in) :: natom,trotter
1365 : real(dp),intent(out) :: eharm,eharm_virial,epot
1366 : !arrays
1367 : real(dp),intent(in) :: etotal_img(trotter),forces(3,natom,trotter)
1368 : real(dp),intent(in) :: xcart(3,natom,trotter)
1369 : real(dp),intent(in) :: spring(natom)
1370 : !Local variables-------------------------------
1371 : !scalars
1372 : integer :: iatom,ii,iimage,iimagep
1373 : !arrays
1374 45 : real(dp),allocatable :: centroid(:,:)
1375 :
1376 : !************************************************************************
1377 :
1378 : !Compute the centroid
1379 135 : ABI_MALLOC(centroid,(3,natom))
1380 645 : centroid=zero
1381 260 : do iimage=1,trotter
1382 830 : do iatom=1,natom
1383 2495 : do ii=1,3
1384 2280 : centroid(ii,iatom)=centroid(ii,iatom)+xcart(ii,iatom,iimage)
1385 : end do
1386 : end do
1387 : end do
1388 645 : centroid=centroid/dble(trotter)
1389 :
1390 : !Potential energy
1391 260 : epot=sum(etotal_img(1:trotter))/dble(trotter)
1392 :
1393 : !Harmonic energy
1394 45 : eharm=zero
1395 260 : do iimage=1,trotter
1396 215 : iimagep=iimage+1;if(iimage==trotter)iimagep=1
1397 830 : do iatom=1,natom
1398 2495 : do ii=1,3
1399 2280 : eharm=eharm+half*spring(iatom)*((xcart(ii,iatom,iimagep)-xcart(ii,iatom,iimage))**2)
1400 : end do
1401 : end do
1402 : end do
1403 :
1404 : !Harmonic energy from virial estimator
1405 45 : eharm_virial=zero
1406 260 : do iimage=1,trotter
1407 830 : do iatom=1,natom
1408 2495 : do ii=1,3
1409 : eharm_virial=eharm_virial-(xcart(ii,iatom,iimage)-centroid(ii,iatom)) &
1410 2280 : & *forces(ii,iatom,iimage)
1411 : end do
1412 : end do
1413 : end do
1414 45 : eharm_virial=eharm_virial/dble(two*trotter)
1415 :
1416 45 : ABI_FREE(centroid)
1417 :
1418 45 : end subroutine pimd_energies
1419 : !!***
1420 :
1421 : !----------------------------------------------------------------------
1422 :
1423 : !!****f* m_pimd/pimd_forces
1424 : !! NAME
1425 : !! pimd_forces
1426 : !!
1427 : !! FUNCTION
1428 : !! Modify forces in order to take into account PIMD contribution
1429 : !!
1430 : !! INPUTS
1431 : !! natom=number of atoms
1432 : !! spring(natom,spring_dim)=spring constants (spring_dim=1 or trotter)
1433 : !! transform=coordinate transformation:
1434 : !! 0: no tranformation
1435 : !! 1: normal mode transformation
1436 : !! 2: staging transformation
1437 : !! trotter=Trotter number
1438 : !! xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell
1439 : !!
1440 : !! OUTPUT
1441 : !!
1442 : !! SIDE EFFECTS
1443 : !! forces(3,natom,trotter)=
1444 : !! at input: forces from electronic calculation
1445 : !! at output: forces from electronic calculation + quantum spring contribution
1446 : !!
1447 : !! SOURCE
1448 :
1449 45 : subroutine pimd_forces(forces,natom,spring,transform,trotter,xcart)
1450 :
1451 : !Arguments ------------------------------------
1452 : !scalars
1453 : integer,intent(in) :: natom,transform,trotter
1454 : !arrays
1455 : real(dp),intent(in) :: xcart(3,natom,trotter)
1456 : real(dp),intent(in) :: spring(:,:)
1457 : real(dp),intent(inout) :: forces(3,natom,trotter)
1458 : !Local variables-------------------------------
1459 : !scalars
1460 : integer :: iatom,ii,iimage,iimagem,iimagep,ispring,natom_spring,nspring
1461 : character(len=500) :: msg
1462 : !arrays
1463 :
1464 : !************************************************************************
1465 :
1466 45 : natom_spring=size(spring,1);nspring=size(spring,2)
1467 45 : if (natom/=natom_spring.or.(nspring/=1.and.nspring/=trotter)) then
1468 0 : msg='Wrong dimensions for array spring !'
1469 0 : ABI_BUG(msg)
1470 : end if
1471 :
1472 45 : if (transform==0) then
1473 120 : do iimage=1,trotter
1474 95 : ispring=min(nspring,iimage)
1475 95 : iimagep=iimage+1; iimagem=iimage-1
1476 95 : if(iimage==trotter) iimagep=1
1477 95 : if(iimage==1) iimagem=trotter
1478 450 : do iatom=1,natom
1479 1415 : do ii=1,3
1480 : forces(ii,iatom,iimage)= &
1481 : & forces(ii,iatom,iimage)/dble(trotter) &
1482 : & - spring(iatom,ispring)*(two*xcart(ii,iatom,iimage)-xcart(ii,iatom,iimagem) &
1483 1320 : & -xcart(ii,iatom,iimagep))
1484 : end do
1485 : end do
1486 : end do
1487 :
1488 : else
1489 140 : do iimage=1,trotter
1490 120 : ispring=min(nspring,iimage)
1491 380 : do iatom=1,natom
1492 1080 : do ii=1,3
1493 : forces(ii,iatom,iimage)= &
1494 : & forces(ii,iatom,iimage)/dble(trotter) &
1495 960 : & - spring(iatom,ispring)*xcart(ii,iatom,iimage)
1496 : end do
1497 : end do
1498 : end do
1499 :
1500 : end if
1501 :
1502 45 : end subroutine pimd_forces
1503 : !!***
1504 :
1505 : !----------------------------------------------------------------------
1506 :
1507 : !!****f* m_pimd/pimd_langevin_forces
1508 : !! NAME
1509 : !! pimd_langevin_forces
1510 : !!
1511 : !! FUNCTION
1512 : !! Compute Langevin contribution to PIMD forces
1513 : !!
1514 : !! INPUTS
1515 : !! alea(3,natom,trotter)=set of random numbers
1516 : !! forces(3,natom,trotter)=forces without Langevin contribution
1517 : !! friction=friction factor
1518 : !! langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter)
1519 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
1520 : !! natom=number of atoms
1521 : !! trotter=Trotter number
1522 : !! vel(3,natom,trotter)=velocities of atoms in each cell
1523 : !!
1524 : !! OUTPUT
1525 : !! forces_langevin(3,natom,trotter)=forces including Langevin contribution
1526 : !!
1527 : !! SIDE EFFECTS
1528 : !!
1529 : !! SOURCE
1530 :
1531 408 : subroutine pimd_langevin_forces(alea,forces,forces_langevin,friction,&
1532 136 : & langev,mass,natom,trotter,vel)
1533 :
1534 : !Arguments ------------------------------------
1535 : !scalars
1536 : integer,intent(in) :: natom,trotter
1537 : real(dp),intent(in) :: friction
1538 : !arrays
1539 : real(dp),intent(in) :: alea(3,natom,trotter),forces(3,natom,trotter)
1540 : real(dp),intent(in) :: langev(:,:),mass(:,:),vel(3,natom,trotter)
1541 : real(dp),intent(out) :: forces_langevin(3,natom,trotter)
1542 : !Local variables-------------------------------
1543 : !scalars
1544 : integer :: iatom,ii,iimage,imass,natom_mass,nmass
1545 : character(len=500) :: msg
1546 : !arrays
1547 :
1548 : !************************************************************************
1549 :
1550 136 : natom_mass=size(mass,1);nmass=size(mass,2)
1551 136 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1552 0 : msg='Wrong dimensions for array mass !'
1553 0 : ABI_BUG(msg)
1554 : end if
1555 :
1556 691 : do iimage=1,trotter
1557 555 : imass=min(nmass,iimage)
1558 2509 : do iatom=1,natom
1559 7827 : do ii=1,3
1560 : forces_langevin(ii,iatom,iimage)=forces(ii,iatom,iimage) &
1561 : & + langev(iatom,imass)*alea(ii,iatom,iimage) &
1562 7272 : & - friction*mass(iatom,imass)*vel(ii,iatom,iimage)
1563 : end do
1564 : end do
1565 : end do
1566 :
1567 136 : end subroutine pimd_langevin_forces
1568 : !!***
1569 :
1570 : !----------------------------------------------------------------------
1571 :
1572 : !!****f* m_pimd/pimd_noseehoover_forces
1573 : !! NAME
1574 : !! pimd_nosehoover_forces
1575 : !!
1576 : !! FUNCTION
1577 : !! Compute Nose-Hoover contribution to PIMD forces
1578 : !! by adding friction force of thermostat number one
1579 : !!
1580 : !! INPUTS
1581 : !! dzeta(3,natom,trotter,nnos)=variables of thermostats, in (atomic time unit)^(-1)
1582 : !! used only when a coordinate transformation is applied (transfom/=0)
1583 : !! forces(3,natom,trotter)=forces without Nose-Hoover contribution
1584 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
1585 : !! natom=number of atoms
1586 : !! nnos=number of thermostats
1587 : !! trotter=Trotter number
1588 : !! vel(3,natom,trotter)=velocities of atoms in each cell
1589 : !!
1590 : !! OUTPUT
1591 : !! forces_nosehoover(3,natom,trotter)=forces including thermostat contribution
1592 : !!
1593 : !! SIDE EFFECTS
1594 : !!
1595 : !! SOURCE
1596 :
1597 132 : subroutine pimd_nosehoover_forces(dzeta,forces,forces_nosehoover,mass,natom,&
1598 132 : & nnos,trotter,vel)
1599 :
1600 : !Arguments ------------------------------------
1601 : !scalars
1602 : integer,intent(in) :: natom,nnos,trotter
1603 : !arrays
1604 : real(dp),intent(in) :: dzeta(3,natom,trotter,nnos),forces(3,natom,trotter)
1605 : real(dp),intent(in) :: vel(3,natom,trotter)
1606 : real(dp),intent(in) :: mass(:,:)
1607 : real(dp),intent(out) :: forces_nosehoover(3,natom,trotter)
1608 : !Local variables-------------------------------
1609 : !scalars
1610 : integer :: iatom,ii,iimage,imass,natom_mass,nmass
1611 : character(len=500) :: msg
1612 : !arrays
1613 :
1614 : !************************************************************************
1615 :
1616 132 : natom_mass=size(mass,1);nmass=size(mass,2)
1617 132 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1618 0 : msg='Wrong dimensions for array mass !'
1619 0 : ABI_BUG(msg)
1620 : end if
1621 :
1622 924 : do iimage=1,trotter
1623 792 : imass=min(nmass,iimage)
1624 2508 : do iatom=1,natom
1625 7128 : do ii=1,3
1626 : forces_nosehoover(ii,iatom,iimage)=forces(ii,iatom,iimage) &
1627 6336 : & - mass(iatom,imass)*dzeta(ii,iatom,iimage,1)*vel(ii,iatom,iimage)
1628 : end do
1629 : end do
1630 : end do
1631 :
1632 132 : end subroutine pimd_nosehoover_forces
1633 : !!***
1634 :
1635 : !----------------------------------------------------------------------
1636 :
1637 : !!****f* m_pimd/pimd_stresses
1638 : !! NAME
1639 : !! pimd_stresses
1640 : !!
1641 : !! FUNCTION
1642 : !! In the case od PIMD, compute the pressure tensor from virial theorem
1643 : !!
1644 : !! INPUTS
1645 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
1646 : !! natom=number of atoms
1647 : !! quantummass(natom)=quantum masses of atoms
1648 : !! stressin(3,3,trotter)=electronic stress tensor for each image
1649 : !! temperature=temperature (could be instantaneous temp. or thermostat temp.)
1650 : !! temperature_therm=thermostat temperature
1651 : !! trotter=Trotter number
1652 : !! vel(3,natom,trotter)=velocities of atoms in each cell
1653 : !! volume=volume of each cell (common to all cells)
1654 : !! xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
1655 : !! temperature=thermostat temperature
1656 : !!
1657 : !! OUTPUT
1658 : !! stress_pimd(3,3,3)=stress tensor for PIMD
1659 : !! First dimension (3) corresponds to 3 different pressure estimators
1660 : !!
1661 : !! SIDE EFFECTS
1662 : !!
1663 : !! SOURCE
1664 :
1665 71 : subroutine pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,&
1666 71 : & temperature,temperature_therm,trotter,vel,volume,xcart)
1667 :
1668 : !Arguments ------------------------------------
1669 : !scalars
1670 : integer,intent(in) :: natom,trotter
1671 : real(dp),intent(in) :: temperature,temperature_therm,volume
1672 : !arrays
1673 : real(dp),intent(in) :: vel(3,natom,trotter),xcart(3,natom,trotter)
1674 : real(dp),intent(in) :: mass(:,:),stressin(3,3,trotter),quantummass(natom)
1675 : real(dp),intent(out) :: stress_pimd(3,3,3)
1676 : !Local variables-------------------------------
1677 : !scalars
1678 : integer :: iatom,ii,iimage,imass,jj,natom_mass,nmass,iimagep
1679 : real(dp) :: stress_tmp(3,3,3),omega2,kt,kt_therm
1680 : character(len=500) :: msg
1681 :
1682 : !************************************************************************
1683 :
1684 71 : natom_mass=size(mass,1);nmass=size(mass,2)
1685 71 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1686 0 : msg='Wrong dimensions for array mass !'
1687 0 : ABI_BUG(msg)
1688 : end if
1689 :
1690 71 : kt=temperature*kb_HaK
1691 71 : kt_therm=temperature_therm*kb_HaK
1692 71 : stress_pimd=zero
1693 :
1694 : !I-PRIMITIVE ESTIMATOR
1695 : !1-Kinetic part
1696 284 : do ii=1,3
1697 284 : stress_pimd(1,ii,ii)=dble(natom)*dble(trotter)*kt/volume
1698 : end do
1699 :
1700 : !2-Potential part
1701 312 : do iimage=1,trotter
1702 1035 : do ii=1,3
1703 3133 : do jj=1,3
1704 2892 : stress_pimd(1,ii,jj)=stress_pimd(1,ii,jj)-stressin(ii,jj,iimage)/dble(trotter)
1705 : !minus to convert stress into pressure
1706 : end do
1707 : end do
1708 : end do
1709 :
1710 : !3-Contribution from springs
1711 71 : omega2=dble(trotter)*kt_therm*kt_therm
1712 312 : do iimage=1,trotter
1713 241 : iimagep=iimage+1
1714 241 : if(iimage==trotter) then
1715 71 : iimagep=1
1716 : end if
1717 1038 : do iatom=1,natom
1718 3145 : do ii=1,3
1719 9438 : do jj=1,3
1720 : stress_pimd(1,ii,jj)=stress_pimd(1,ii,jj)-quantummass(iatom)*omega2* &
1721 : & (xcart(ii,iatom,iimagep)-xcart(ii,iatom,iimage))* &
1722 8712 : & (xcart(jj,iatom,iimagep)-xcart(jj,iatom,iimage))/volume
1723 : end do
1724 : end do
1725 : end do
1726 : end do
1727 :
1728 : !II-Average of classical pressures
1729 : !1-Kinetic part
1730 312 : do iimage=1,trotter
1731 241 : imass=min(nmass,iimage)
1732 1038 : do iatom=1,natom
1733 3145 : do ii=1,3
1734 9438 : do jj=1,3
1735 : stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)+ &
1736 8712 : & mass(iatom,imass)*vel(ii,iatom,iimage)*vel(jj,iatom,iimage)/volume
1737 : end do
1738 : end do
1739 : end do
1740 : end do
1741 :
1742 : !2-Contribution from electronic stress
1743 312 : do iimage=1,trotter
1744 1035 : do ii=1,3
1745 3133 : do jj=1,3
1746 2892 : stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)-stressin(ii,jj,iimage)
1747 : end do
1748 : end do
1749 : end do
1750 284 : do ii=1,3
1751 923 : do jj=1,3
1752 852 : stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)/dble(trotter)
1753 : end do
1754 : end do
1755 :
1756 : !III-pressure from VIRIAL estimator: stress_pimd(3,:,:)
1757 : !1-kinetic part
1758 : ! do ii=1,3
1759 : ! stress_pimd(3,ii,ii)=dfloat(natom)*kt/volume
1760 : ! end do
1761 :
1762 : !Symmetrize internal pressure
1763 71 : stress_tmp=stress_pimd
1764 284 : stress_pimd(:,2,1)=half*(stress_tmp(:,1,2)+stress_tmp(:,2,1))
1765 284 : stress_pimd(:,1,2)=half*(stress_tmp(:,1,2)+stress_tmp(:,2,1))
1766 284 : stress_pimd(:,3,1)=half*(stress_tmp(:,1,3)+stress_tmp(:,3,1))
1767 284 : stress_pimd(:,1,3)=half*(stress_tmp(:,1,3)+stress_tmp(:,3,1))
1768 284 : stress_pimd(:,2,3)=half*(stress_tmp(:,3,2)+stress_tmp(:,2,3))
1769 284 : stress_pimd(:,3,2)=half*(stress_tmp(:,3,2)+stress_tmp(:,2,3))
1770 :
1771 71 : end subroutine pimd_stresses
1772 : !!***
1773 :
1774 : !----------------------------------------------------------------------
1775 :
1776 : !!****f* m_pimd/pimd_diff_stress
1777 : !! NAME
1778 : !! pimd_diff_stress
1779 : !!
1780 : !! FUNCTION
1781 : !! Compute the difference between the stress tensor and the stress target
1782 : !!
1783 : !! INPUTS
1784 : !! stress_pimd(3,3,3)=stress tensor for PIMD
1785 : !! Last dimension (3) corresponds to 3 different pressure estimators
1786 : !! stress_target(6)=stress target
1787 : !!
1788 : !! OUTPUT
1789 : !! pimd_diff_stress(3,3,3)=difference between stresses and stress target
1790 : !! First dimension (3) corresponds to 3 different pressure estimators
1791 : !!
1792 : !! SIDE EFFECTS
1793 : !!
1794 : !! SOURCE
1795 :
1796 26 : function pimd_diff_stress(stress_pimd,stress_target)
1797 :
1798 : !Arguments ------------------------------------
1799 : !scalars
1800 : !arrays
1801 : real(dp),intent(in) :: stress_pimd(3,3,3),stress_target(6)
1802 : real(dp) :: pimd_diff_stress(3,3)
1803 : !Local variables-------------------------------
1804 : !scalars
1805 : !arrays
1806 : real(dp) :: stress_pimd2(3,3)
1807 :
1808 : !************************************************************************
1809 :
1810 : !Choice: the primitive estimator for pressure is chosen
1811 338 : stress_pimd2(:,:) = stress_pimd(1, :, :)
1812 :
1813 : !+stress_target instead of - because it is translated from stress to pressure tensor
1814 26 : pimd_diff_stress(1,1)=stress_pimd2(1,1)+stress_target(1)
1815 26 : pimd_diff_stress(2,2)=stress_pimd2(2,2)+stress_target(2)
1816 26 : pimd_diff_stress(3,3)=stress_pimd2(3,3)+stress_target(3)
1817 26 : pimd_diff_stress(2,3)=stress_pimd2(2,3)+stress_target(4)
1818 26 : pimd_diff_stress(3,2)=stress_pimd2(3,2)+stress_target(4)
1819 26 : pimd_diff_stress(1,3)=stress_pimd2(1,3)+stress_target(5)
1820 26 : pimd_diff_stress(3,1)=stress_pimd2(3,1)+stress_target(5)
1821 26 : pimd_diff_stress(1,2)=stress_pimd2(1,2)+stress_target(6)
1822 26 : pimd_diff_stress(2,1)=stress_pimd2(2,1)+stress_target(6)
1823 :
1824 : end function pimd_diff_stress
1825 : !!***
1826 :
1827 : !----------------------------------------------------------------------
1828 :
1829 : !!****f* m_pimd/pimd_predict_taylor
1830 : !! NAME
1831 : !! pimd_predict_taylor
1832 : !!
1833 : !! FUNCTION
1834 : !! Predict new atomic positions using a Taylor algorithm (first time step) - for PIMD
1835 : !!
1836 : !! INPUTS
1837 : !! dtion=time step
1838 : !! forces(3,natom,trotter)=PIMD forces on atoms in each cell
1839 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
1840 : !! natom=number of atoms
1841 : !! trotter=Trotter number
1842 : !! vel(3,natom,trotter)=velocities of atoms in each cell
1843 : !! xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
1844 : !!
1845 : !! OUTPUT
1846 : !! xcart_next(3,natom,trotter)=cartesian coordinates of atoms in each cell at t+dt
1847 : !!
1848 : !! SIDE EFFECTS
1849 : !!
1850 : !! SOURCE
1851 :
1852 7 : subroutine pimd_predict_taylor(dtion,forces,mass,natom,trotter,vel,xcart,xcart_next)
1853 :
1854 : !Arguments ------------------------------------
1855 : !scalars
1856 : integer,intent(in) :: natom,trotter
1857 : real(dp),intent(in) :: dtion
1858 : !arrays
1859 : real(dp),intent(in) :: forces(3,natom,trotter),vel(3,natom,trotter),xcart(3,natom,trotter)
1860 : real(dp),intent(in) :: mass(:,:)
1861 : real(dp),intent(out) :: xcart_next(3,natom,trotter)
1862 : !Local variables-------------------------------
1863 : !scalars
1864 : integer :: iatom,ii,iimage,imass,natom_mass,nmass
1865 : character(len=500) :: msg
1866 : !arrays
1867 :
1868 : !************************************************************************
1869 :
1870 7 : natom_mass=size(mass,1);nmass=size(mass,2)
1871 7 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1872 0 : msg='Wrong dimensions for array mass !'
1873 0 : ABI_BUG(msg)
1874 : end if
1875 :
1876 44 : do iimage=1,trotter
1877 37 : imass=min(nmass,iimage)
1878 122 : do iatom=1,natom
1879 349 : do ii=1,3
1880 : xcart_next(ii,iatom,iimage)=xcart(ii,iatom,iimage) &
1881 : & + half*dtion*dtion*forces(ii,iatom,iimage)/mass(iatom,imass) &
1882 312 : & + dtion*vel(ii,iatom,iimage)
1883 : end do
1884 : end do
1885 : end do
1886 :
1887 7 : end subroutine pimd_predict_taylor
1888 : !!***
1889 :
1890 : !----------------------------------------------------------------------
1891 :
1892 : !!****f* m_pimd/pimd_predict_verlet
1893 : !! NAME
1894 : !! pimd_predict_verlet
1895 : !!
1896 : !! FUNCTION
1897 : !! Predict new atomic positions using a Verlet algorithm - for PIMD
1898 : !!
1899 : !! INPUTS
1900 : !! dtion=time step
1901 : !! forces(3,natom,trotter)=PIMD forces on atoms in each cell
1902 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
1903 : !! natom=number of atoms
1904 : !! trotter=Trotter number
1905 : !! xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
1906 : !! xcart_prev(3,natom,trotter)=cartesian coordinates of atoms in each cell at t-dt
1907 : !!
1908 : !! OUTPUT
1909 : !! xcart_next(3,natom,trotter)=cartesian coordinates of atoms in each cell at t+dt
1910 : !!
1911 : !! SIDE EFFECTS
1912 : !!
1913 : !! SOURCE
1914 :
1915 235 : subroutine pimd_predict_verlet(dtion,forces,mass,natom,trotter,xcart,xcart_next,xcart_prev)
1916 :
1917 : !Arguments ------------------------------------
1918 : !scalars
1919 : integer,intent(in) :: natom,trotter
1920 : real(dp),intent(in) :: dtion
1921 : !arrays
1922 : real(dp),intent(in) :: forces(3,natom,trotter)
1923 : real(dp),intent(in) :: xcart(3,natom,trotter),xcart_prev(3,natom,trotter)
1924 : real(dp),intent(in) :: mass(:,:)
1925 : real(dp),intent(out) :: xcart_next(3,natom,trotter)
1926 : !Local variables-------------------------------
1927 : !scalars
1928 : integer :: iatom,ii,iimage,imass,natom_mass,nmass
1929 : character(len=500) :: msg
1930 : !arrays
1931 :
1932 : !************************************************************************
1933 :
1934 235 : natom_mass=size(mass,1);nmass=size(mass,2)
1935 235 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1936 0 : msg='Wrong dimensions for array mass !'
1937 0 : ABI_BUG(msg)
1938 : end if
1939 :
1940 1519 : do iimage=1,trotter
1941 1284 : imass=min(nmass,iimage)
1942 4687 : do iatom=1,natom
1943 13956 : do ii=1,3
1944 : xcart_next(ii,iatom,iimage)= &
1945 : & two*xcart(ii,iatom,iimage) &
1946 : & - xcart_prev(ii,iatom,iimage) &
1947 12672 : & + dtion*dtion*forces(ii,iatom,iimage)/mass(iatom,imass)
1948 : end do
1949 : end do
1950 : end do
1951 :
1952 235 : end subroutine pimd_predict_verlet
1953 : !!***
1954 :
1955 : !----------------------------------------------------------------------
1956 :
1957 : !!****f* m_pimd/pimd_predict_vel
1958 : !! NAME
1959 : !! pimd_predict_vel
1960 : !!
1961 : !! FUNCTION
1962 : !! Predict an estimation of the velocities at next time time step
1963 : !! from the positions of 3 consecutive time steps
1964 : !!
1965 : !! INPUTS
1966 : !! dtion=time step
1967 : !! itime=time step index
1968 : !! natom=number of atoms
1969 : !! trotter=Trotter number
1970 : !! pitransform=integer selecting the transformation of coordinates
1971 : !! xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
1972 : !! xcart_next(3,natom,trotter)=cartesian coordinates of atoms in each cell at t+dt
1973 : !! xcart_prev(3,natom,trotter)=cartesian coordinates of atoms in each cell at t-dt
1974 : !!
1975 : !! OUTPUT
1976 : !! vel(3,natom,trotter)=estimation of velocities at t+dt
1977 : !!
1978 : !! SIDE EFFECTS
1979 : !!
1980 : !! SOURCE
1981 :
1982 44 : subroutine pimd_predict_vel(dtion,itime,natom,trotter,pitransform,xcart,xcart_next,xcart_prev,vel)
1983 :
1984 : !Arguments ------------------------------------
1985 : !scalars
1986 : integer,intent(in) :: itime,natom,pitransform,trotter
1987 : real(dp),intent(in) :: dtion
1988 : !arrays
1989 : real(dp),intent(inout) :: xcart(3,natom,trotter),xcart_next(3,natom,trotter),xcart_prev(3,natom,trotter)
1990 : real(dp),intent(out) :: vel(3,natom,trotter)
1991 : !Local variables-------------------------------
1992 : !scalars
1993 : integer :: iatom,ii,iimage
1994 :
1995 : !************************************************************************
1996 :
1997 : !Set the coordinates in the transformed representation
1998 44 : if (pitransform>0) then
1999 10 : call pimd_coord_transform(xcart_next,1,natom,pitransform,trotter)
2000 10 : call pimd_coord_transform(xcart,1,natom,pitransform,trotter)
2001 10 : call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter)
2002 : end if
2003 :
2004 44 : if (itime>1) then
2005 214 : do iimage=1,trotter
2006 700 : do iatom=1,natom
2007 2121 : do ii=1,3
2008 : vel(ii,iatom,iimage)= &
2009 : & (three*xcart_next(ii,iatom,iimage) &
2010 : & -four*xcart(ii,iatom,iimage) &
2011 : & +xcart_prev(ii,iatom,iimage)) &
2012 1944 : & / (two*dtion)
2013 : end do
2014 : end do
2015 : end do
2016 : else
2017 44 : do iimage=1,trotter
2018 122 : do iatom=1,natom
2019 349 : do ii=1,3
2020 312 : vel(ii,iatom,iimage)=(xcart_next(ii,iatom,iimage)-xcart(ii,iatom,iimage))/dtion
2021 : end do
2022 : end do
2023 : end do
2024 : end if
2025 :
2026 : !Set back the coordinates
2027 44 : if (pitransform>0) then
2028 10 : call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter)
2029 10 : call pimd_coord_transform(xcart,-1,natom,pitransform,trotter)
2030 10 : call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter)
2031 : end if
2032 :
2033 44 : end subroutine pimd_predict_vel
2034 : !!***
2035 :
2036 : !----------------------------------------------------------------------
2037 :
2038 : !!****f* m_pimd/pimd_predict_velcell
2039 : !! NAME
2040 : !! pimd_predict_velcell
2041 : !!
2042 : !! FUNCTION
2043 : !! Predict an estimation of the cell velocities at next time time step
2044 : !! from the cell at 3 consecutive time steps
2045 : !!
2046 : !! INPUTS
2047 : !! dtion=time step
2048 : !! itime=time step index
2049 : !! pitransform=integer selecting the transformation of coordinates
2050 : !! rprimd(3,3)=cell coordinates at t
2051 : !! rprimd_next(3,3)=cell coordinates at t+dt
2052 : !! rprimd_prev(3,3)=cell coordinates at t-dt
2053 : !!
2054 : !! OUTPUT
2055 : !! vel_cell(3,3)=estimation of cell velocities at t+dt
2056 : !!
2057 : !! SIDE EFFECTS
2058 : !!
2059 : !! SOURCE
2060 :
2061 4 : subroutine pimd_predict_velcell(dtion,itime,rprimd,rprimd_next,rprimd_prev,vel_cell)
2062 :
2063 : !Arguments ------------------------------------
2064 : !scalars
2065 : integer,intent(in) :: itime
2066 : real(dp),intent(in) :: dtion
2067 : !arrays
2068 : real(dp),intent(in) :: rprimd(3,3),rprimd_next(3,3),rprimd_prev(3,3)
2069 : real(dp),intent(out) :: vel_cell(3,3)
2070 : !Local variables-------------------------------
2071 : !scalars
2072 : integer :: ii,jj
2073 :
2074 : !************************************************************************
2075 :
2076 4 : if (itime>1) then
2077 16 : do jj=1,3
2078 52 : do ii=1,3
2079 : vel_cell(ii,jj)= &
2080 48 : & (three*rprimd_next(ii,jj)-four*rprimd(ii,jj)+rprimd_prev(ii,jj))/(two*dtion)
2081 : end do
2082 : end do
2083 : else
2084 0 : do jj=1,3
2085 0 : do ii=1,3
2086 0 : vel_cell(ii,jj)=(rprimd_next(ii,jj)-rprimd(ii,jj))/dtion
2087 : end do
2088 : end do
2089 : end if
2090 :
2091 4 : end subroutine pimd_predict_velcell
2092 : !!***
2093 :
2094 : !----------------------------------------------------------------------
2095 :
2096 : !!****f* m_pimd/pimd_nosehoover_propagate
2097 : !! NAME
2098 : !! pimd_nosehoover_propagate
2099 : !!
2100 : !! FUNCTION
2101 : !! Propagate thermostat variables (Nose-Hoover algorithm) - for PIMD
2102 : !!
2103 : !! INPUTS
2104 : !! dtion=time step
2105 : !! dzeta(3,natom,trotter,nnos)=time derivative of zeta at time t
2106 : !! itimimage=index of time step
2107 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
2108 : !! natom=number of atoms
2109 : !! nnos=number of thermostats
2110 : !! qmass(nnos)=masses of thermostats
2111 : !! temperature=temperature
2112 : !! transform=coordinate transformation:
2113 : !! 0: no tranformation
2114 : !! 1: normal mode transformation
2115 : !! 2: staging transformation
2116 : !! trotter=Trotter number
2117 : !! vel(3,natom,trotter)=velocities of atoms in each cell
2118 : !! zeta(3,natom,trotter,nnos)=variables of thermostats, in (atomic time unit)^(-1)
2119 : !! used only when no coordinate transformation is applied (transfom==0)
2120 : !! at time t
2121 : !! zeta_prev(3,natom,trotter,nnos)=previous value of zeta (t-dt)
2122 : !!
2123 : !! OUTPUT
2124 : !! zeta_next(3,natom,trotter,nnos)=next value of zeta (t+dt)
2125 : !!
2126 : !! SIDE EFFECTS
2127 : !!
2128 : !! SOURCE
2129 :
2130 132 : subroutine pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,temperature,&
2131 132 : & trotter,vel,zeta,zeta_next,zeta_prev,itimimage,transform)
2132 :
2133 : !Arguments ------------------------------------
2134 : !scalars
2135 : integer,intent(in) :: itimimage,natom,nnos,transform,trotter
2136 : real(dp),intent(in) :: dtion,temperature
2137 : !arrays
2138 : real(dp),intent(in) :: qmass(nnos),vel(3,natom,trotter)
2139 : real(dp),intent(in) :: mass(:,:)
2140 : real(dp),intent(in) :: dzeta(3,natom,trotter,nnos),zeta_prev(3,natom,trotter,nnos)
2141 : real(dp),intent(in) :: zeta(3,natom,trotter,nnos)
2142 : real(dp),intent(out) :: zeta_next(3,natom,trotter,nnos)
2143 : !Local variables-------------------------------
2144 : !scalars
2145 : integer :: iatom,ii,iimage,inos,natom_mass,nmass
2146 : real(dp) :: kt
2147 : character(len=500) :: msg
2148 : !arrays
2149 132 : real(dp),allocatable :: thermforces(:,:,:,:)
2150 :
2151 : !************************************************************************
2152 :
2153 132 : natom_mass=size(mass,1);nmass=size(mass,2)
2154 132 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
2155 0 : msg='Wrong dimensions for array mass !'
2156 0 : ABI_BUG(msg)
2157 : end if
2158 132 : if (nnos<3) then
2159 0 : msg='Not available for nnos<3 !'
2160 0 : ABI_BUG(msg)
2161 : end if
2162 :
2163 132 : kt=temperature*kb_HaK
2164 660 : ABI_MALLOC(thermforces,(3,natom,trotter,nnos))
2165 :
2166 : !forces on the thermostats
2167 792 : do inos=1,nnos
2168 792 : if(inos==1)then
2169 924 : do iimage=1,trotter
2170 2508 : do iatom=1,natom
2171 7128 : do ii=1,3
2172 1584 : select case(transform)
2173 : case(0)
2174 1512 : thermforces(ii,iatom,iimage,1)=mass(iatom,1)*(vel(ii,iatom,iimage)**2)-kt
2175 : case(1)
2176 1620 : thermforces(ii,iatom,iimage,1)=mass(iatom,iimage)*(vel(ii,iatom,iimage)**2)-kt
2177 : case(2)
2178 4752 : thermforces(ii,iatom,iimage,1)=mass(iatom,iimage)*(vel(ii,iatom,iimage)**2)-kt
2179 : end select
2180 : end do
2181 : end do
2182 : end do
2183 528 : elseif(inos==(nnos-1))then
2184 924 : do iimage=1,trotter
2185 2508 : do iatom=1,natom
2186 7128 : do ii=1,3
2187 : thermforces(ii,iatom,iimage,inos)=&
2188 : & (qmass(nnos-2)*(dzeta(ii,iatom,iimage,nnos-2)**2)-kt) + &
2189 6336 : & (qmass(nnos )*(dzeta(ii,iatom,iimage,nnos) **2)-kt)
2190 : end do
2191 : end do
2192 : end do
2193 396 : elseif(inos==nnos)then
2194 924 : do iimage=1,trotter
2195 2508 : do iatom=1,natom
2196 7128 : do ii=1,3
2197 : thermforces(ii,iatom,iimage,inos)=&
2198 6336 : & qmass(nnos-1)*(dzeta(ii,iatom,iimage,nnos-1)**2)-kt
2199 : end do
2200 : end do
2201 : end do
2202 : else
2203 1848 : do iimage=1,trotter
2204 5016 : do iatom=1,natom
2205 14256 : do ii=1,3
2206 : thermforces(ii,iatom,iimage,inos)=&
2207 12672 : & qmass(inos-1)*(dzeta(ii,iatom,iimage,inos-1)**2)-kt
2208 : end do
2209 : end do
2210 : end do
2211 : end if
2212 : end do
2213 :
2214 132 : select case(itimimage)
2215 : case(1) !taylor
2216 :
2217 18 : do inos=1,nnos
2218 18 : if(inos==1)then
2219 : zeta_next(:,:,:,1)=zeta(:,:,:,1)+ dzeta(:,:,:,1)*dtion + &
2220 : & (thermforces(:,:,:,1)-qmass(1)*dzeta(:,:,:,1)*dzeta(:,:,:,2))* &
2221 165 : & dtion*dtion/(two*qmass(1))
2222 12 : elseif(inos==(nnos-1))then
2223 : zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
2224 : & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,nnos))* &
2225 165 : & dtion*dtion/(two*qmass(inos))
2226 9 : elseif(inos==nnos)then
2227 : zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
2228 : & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,nnos-1)*dzeta(:,:,:,nnos))* &
2229 165 : & dtion*dtion/(two*qmass(inos))
2230 : else
2231 : zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
2232 : & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,inos+1))* &
2233 330 : & dtion*dtion/(two*qmass(inos))
2234 : end if
2235 : end do
2236 :
2237 : case default !verlet
2238 :
2239 906 : do inos=1,nnos
2240 774 : if(inos==1)then
2241 : zeta_next(:,:,:,1)=two*zeta(:,:,:,1) - zeta_prev(:,:,:,1) + &
2242 : & (thermforces(:,:,:,1)-qmass(1)*dzeta(:,:,:,1)*dzeta(:,:,:,2))* &
2243 7095 : & dtion*dtion/qmass(1)
2244 516 : elseif(inos==(nnos-1))then
2245 : zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2246 : & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,nnos))* &
2247 7095 : & dtion*dtion/qmass(inos)
2248 387 : elseif(inos==nnos)then
2249 : zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2250 : & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,nnos-1)*dzeta(:,:,:,nnos))* &
2251 7095 : & dtion*dtion/qmass(inos)
2252 : else
2253 : zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2254 : & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,inos+1))* &
2255 14190 : & dtion*dtion/qmass(inos)
2256 : end if
2257 : end do
2258 :
2259 : end select
2260 :
2261 132 : ABI_FREE(thermforces)
2262 :
2263 132 : end subroutine pimd_nosehoover_propagate
2264 : !!***
2265 :
2266 : !----------------------------------------------------------------------
2267 :
2268 : !!****f* m_pimd/pimd_coord_transform
2269 : !! NAME
2270 : !! pimd_coord_transform
2271 : !!
2272 : !! FUNCTION
2273 : !! Apply a coordinate transformation on a given vector field
2274 : !! (defined for each atom in in cell) - for PIMD
2275 : !! Possible choices for the transformation:
2276 : !! 0: no transformation
2277 : !! 1: normal mode tranformation
2278 : !! 2: staging transformation
2279 : !!
2280 : !! INPUTS
2281 : !! ioption=option given the direction of the transformation
2282 : !! +1: from primitive coordinates to transformed coordinates
2283 : !! -1: from transformed coordinates to primitive coordinates
2284 : !! natom=number of atoms
2285 : !! transform=coordinate transformation:
2286 : !! 0: no tranformation
2287 : !! 1: normal mode transformation
2288 : !! 2: staging transformation
2289 : !! trotter=Trotter number
2290 : !!
2291 : !! OUTPUT
2292 : !!
2293 : !! SIDE EFFECTS
2294 : !! array(3,natom,trotter)=array to be transformed
2295 : !!
2296 : !! SOURCE
2297 :
2298 253 : subroutine pimd_coord_transform(array,ioption,natom,transform,trotter)
2299 :
2300 : !Arguments ------------------------------------
2301 : !scalars
2302 : integer,intent(in) :: ioption,natom,transform,trotter
2303 : !arrays
2304 : real(dp),intent(inout) :: array(3,natom,trotter)
2305 : !Local variables-------------------------------
2306 : !scalars
2307 : integer :: iatom,ii,iimage,iimagem,iimagep,iimage2,ll
2308 : !arrays
2309 253 : real(dp),allocatable :: array_temp(:,:,:)
2310 253 : real(dp),allocatable :: nrm(:,:)
2311 :
2312 : !************************************************************************
2313 :
2314 : !=== No transformation ===================================================
2315 253 : if (transform==0) then
2316 :
2317 : return
2318 :
2319 : !=== Normal mode transformation ==========================================
2320 156 : else if (transform==1) then
2321 :
2322 : ! ---------------- From primitive to transformed coordinates ------------
2323 24 : if (ioption==+1) then
2324 36 : ABI_MALLOC(array_temp,(3,natom,trotter))
2325 :
2326 81 : array_temp(:,:,1)=zero
2327 63 : do iimage=1,trotter
2328 495 : array_temp(:,:,1)=array_temp(:,:,1)+array(:,:,iimage)
2329 : end do
2330 :
2331 27 : do ll=2,trotter/2
2332 162 : array_temp(:,:,2*ll-2)=zero
2333 162 : array_temp(:,:,2*ll-1)=zero
2334 135 : do iimage=1,trotter
2335 : array_temp(:,:,2*ll-2)=array_temp(:,:,2*ll-2)+array(:,:,iimage)* &
2336 972 : & cos(two_pi*dble((iimage-1)*(ll-1))/dble(trotter))
2337 : array_temp(:,:,2*ll-1)=array_temp(:,:,2*ll-1)-array(:,:,iimage)* &
2338 990 : & sin(two_pi*dble((iimage-1)*(ll-1))/dble(trotter))
2339 : end do
2340 : end do
2341 :
2342 81 : array_temp(:,:,trotter)=zero
2343 63 : do iimage=1,trotter
2344 : array_temp(:,:,trotter)=array_temp(:,:,trotter)+&
2345 495 : & array(:,:,iimage)*dble((-1)**(iimage-1))
2346 : end do
2347 :
2348 495 : array=array_temp/dble(trotter)
2349 :
2350 9 : ABI_FREE(array_temp)
2351 :
2352 : ! ---------------- From transformed to primitive coordinates ------------
2353 15 : else if (ioption==-1) then
2354 :
2355 60 : ABI_MALLOC(array_temp,(3,natom,trotter)) !real part
2356 60 : ABI_MALLOC(nrm,(trotter,trotter)) !real part
2357 :
2358 105 : do iimage=1,trotter
2359 90 : nrm(iimage,1)=one
2360 105 : nrm(iimage,trotter)=dble((-1)**(iimage-1))
2361 : end do
2362 :
2363 45 : do ll=2,trotter/2
2364 225 : do iimage=1,trotter
2365 : nrm(iimage,2*ll-2)= two*cos(two_pi*dble((ll-1)* &
2366 180 : & (iimage-1))/dble(trotter))
2367 : nrm(iimage,2*ll-1)=-two*sin(two_pi*dble((ll-1)* &
2368 210 : & (iimage-1))/dble(trotter))
2369 : end do
2370 : end do
2371 :
2372 105 : do iimage=1,trotter
2373 810 : array_temp(:,:,iimage)=zero
2374 645 : do ll=1,trotter
2375 4950 : array_temp(:,:,iimage)=array_temp(:,:,iimage)+nrm(iimage,ll)*array(:,:,ll)
2376 : end do
2377 : end do
2378 825 : array=array_temp
2379 :
2380 15 : ABI_FREE(array_temp)
2381 15 : ABI_FREE(nrm)
2382 :
2383 : end if ! ioption
2384 :
2385 : !=== Staging transformation ==============================================
2386 132 : else if (transform==2) then
2387 :
2388 : ! ---------------- From primitive to transformed coordinates ------------
2389 132 : if (ioption==+1) then
2390 :
2391 228 : ABI_MALLOC(array_temp,(3,natom,trotter))
2392 3135 : array_temp=zero
2393 399 : do iimage=1,trotter
2394 342 : iimagep=iimage+1;if(iimage==trotter) iimagep=1
2395 342 : iimagem=iimage-1;if(iimage==1) iimagem=trotter
2396 1083 : do iatom=1,natom
2397 3078 : do ii=1,3
2398 2736 : array_temp(ii,iatom,iimage)=(dble(iimagem)*array(ii,iatom,iimagep)+array(ii,iatom,1))/dble(iimage)
2399 : end do
2400 : end do
2401 : end do
2402 57 : if (trotter>1) then
2403 342 : do iimage=2,trotter
2404 912 : do iatom=1,natom
2405 2565 : do ii=1,3
2406 2280 : array(ii,iatom,iimage)=array(ii,iatom,iimage)-array_temp(ii,iatom,iimage)
2407 : end do
2408 : end do
2409 : end do
2410 : end if
2411 57 : ABI_FREE(array_temp)
2412 :
2413 : ! ---------------- From transformed to primitive coordinates ------------
2414 75 : else if (ioption==-1) then
2415 :
2416 300 : ABI_MALLOC(array_temp,(3,natom,trotter))
2417 4125 : array_temp=zero
2418 525 : do iimage=1,trotter
2419 1425 : do iatom=1,natom
2420 4050 : do ii=1,3
2421 3600 : array_temp(ii,iatom,iimage)=array(ii,iatom,1)
2422 : end do
2423 : end do
2424 : end do
2425 75 : if (trotter>1) then
2426 450 : do iimage=2,trotter
2427 1200 : do iatom=1,natom
2428 3375 : do ii=1,3
2429 9750 : do iimage2=iimage,trotter
2430 : array_temp(ii,iatom,iimage)=array_temp(ii,iatom,iimage) &
2431 9000 : & +array(ii,iatom,iimage2)*(dble(iimage-1))/(dble(iimage2-1))
2432 : end do
2433 : end do
2434 : end do
2435 : end do
2436 : end if
2437 4125 : array(:,:,:)=array_temp(:,:,:)
2438 75 : ABI_FREE(array_temp)
2439 :
2440 : end if ! ioption
2441 :
2442 : end if ! transform
2443 :
2444 : end subroutine pimd_coord_transform
2445 : !!***
2446 :
2447 : !----------------------------------------------------------------------
2448 :
2449 : !!****f* m_pimd/pimd_force_transform
2450 : !! NAME
2451 : !! pimd_force_transform
2452 : !!
2453 : !! FUNCTION
2454 : !! Apply a coordinate transformation on forces (defined for each atom in in cell) - for PIMD
2455 : !! Possible choices for the transformation:
2456 : !! 0: no transformation
2457 : !! 1: normal mode tranformation
2458 : !! 2: staging transformation
2459 : !!
2460 : !! INPUTS
2461 : !! ioption=option given the direction of the transformation
2462 : !! +1: from primitive coordinates to transformed coordinates
2463 : !! -1: from transformed coordinates to primitive coordinates
2464 : !! natom=number of atoms
2465 : !! transform=coordinate transformation:
2466 : !! 0: no tranformation
2467 : !! 1: normal mode transformation
2468 : !! 2: staging transformation
2469 : !! trotter=Trotter number
2470 : !!
2471 : !! OUTPUT
2472 : !!
2473 : !! SIDE EFFECTS
2474 : !! forces(3,natom,trotter)=array containing forces
2475 : !!
2476 : !! NOTES
2477 : !! Back transformation (ioption=-1) not implemented !
2478 : !!
2479 : !! SOURCE
2480 :
2481 40 : subroutine pimd_force_transform(forces,ioption,natom,transform,trotter)
2482 :
2483 : !Arguments ------------------------------------
2484 : !scalars
2485 : integer,intent(in) :: ioption,natom,transform,trotter
2486 : !arrays
2487 : real(dp),intent(inout) :: forces(3,natom,trotter)
2488 : !Local variables-------------------------------
2489 : !scalars
2490 : integer :: iatom,ii,iimage,ll
2491 : character(len=500) :: msg
2492 : !arrays
2493 40 : real(dp),allocatable :: forces_temp(:,:,:),nrm(:,:)
2494 :
2495 : !************************************************************************
2496 :
2497 40 : if (ioption==-1) then
2498 0 : msg='Back transformation not implemented !'
2499 0 : ABI_BUG(msg)
2500 : end if
2501 :
2502 : !=== No transformation ===================================================
2503 40 : if (transform==0) then
2504 :
2505 : return
2506 :
2507 : !=== Normal mode transformation ==========================================
2508 20 : else if (transform==1) then
2509 :
2510 : !normal mode forces
2511 20 : ABI_MALLOC(forces_temp,(3,natom,trotter))
2512 20 : ABI_MALLOC(nrm,(trotter,trotter))
2513 :
2514 35 : do iimage=1,trotter
2515 30 : nrm(iimage,1)=one
2516 35 : nrm(iimage,trotter)=dble((-1)**(iimage-1))
2517 : end do
2518 :
2519 15 : do ll=2,trotter/2
2520 75 : do iimage=1,trotter
2521 : nrm(iimage,2*ll-2)= two*cos(two_pi*dble((ll-1)* &
2522 60 : & (iimage-1))/dble(trotter))
2523 : nrm(iimage,2*ll-1)=-two*sin(two_pi*dble((ll-1)* &
2524 70 : & (iimage-1))/dble(trotter))
2525 : end do
2526 : end do
2527 :
2528 35 : do ll=1,trotter
2529 270 : forces_temp(:,:,ll)=zero
2530 215 : do iimage=1,trotter
2531 1650 : forces_temp(:,:,ll)=forces_temp(:,:,ll)+nrm(iimage,ll)*forces(:,:,iimage)
2532 : end do
2533 : end do
2534 :
2535 275 : forces=forces_temp
2536 :
2537 5 : ABI_FREE(forces_temp)
2538 5 : ABI_FREE(nrm)
2539 :
2540 : !=== Staging transformation ==============================================
2541 15 : else if (transform==2) then
2542 :
2543 : !staging forces
2544 60 : ABI_MALLOC(forces_temp,(3,natom,trotter))
2545 825 : forces_temp=zero
2546 105 : do iimage=1,trotter
2547 285 : do iatom=1,natom
2548 810 : do ii=1,3
2549 720 : forces_temp(ii,iatom,1)=forces_temp(ii,iatom,1)+forces(ii,iatom,iimage)
2550 : end do
2551 : end do
2552 : end do
2553 15 : if (trotter>1) then
2554 90 : do iimage=2,trotter
2555 240 : do iatom=1,natom
2556 675 : do ii=1,3
2557 : forces_temp(ii,iatom,iimage)=forces(ii,iatom,iimage) &
2558 600 : & +forces_temp(ii,iatom,iimage-1)*(dble(iimage-2)/dble(iimage-1))
2559 : end do
2560 : end do
2561 : end do
2562 : end if
2563 825 : forces=forces_temp
2564 15 : ABI_FREE(forces_temp)
2565 :
2566 : end if ! transform
2567 :
2568 : end subroutine pimd_force_transform
2569 : !!***
2570 :
2571 : !----------------------------------------------------------------------
2572 :
2573 : !!****f* m_pimd/pimd_apply_constraint
2574 : !! NAME
2575 : !! pimd_apply_constraint
2576 : !!
2577 : !! FUNCTION
2578 : !! Modify forces to take into account an holonomic constraint
2579 : !! according to "pimd_constraint" parameter
2580 : !! Available constraints:
2581 : !! 0: no constraint
2582 : !! 1: linear combination of coordinates
2583 : !!
2584 : !! INPUTS
2585 : !! constraint=type of constraint to be applied
2586 : !! mass(natom)=fictitious masses of atoms
2587 : !! natom=number of atoms
2588 : !! trotter=Trotter number
2589 : !! wtatcon(3,natom)=weights for atomic constraints
2590 : !! xcart(3,natom,trotter)=cartesian coordinates of atoms
2591 : !!
2592 : !! OUTPUT
2593 : !! constraint_output(2)=several (real) data to be output
2594 : !! when a constraint has been applied
2595 : !!
2596 : !! SIDE EFFECTS
2597 : !! forces(3,natom,trotter)=array containing forces
2598 : !!
2599 : !! SOURCE
2600 :
2601 268 : subroutine pimd_apply_constraint(constraint,constraint_output,forces,mass,natom,&
2602 268 : & trotter,wtatcon,xcart)
2603 :
2604 : !Arguments ------------------------------------
2605 : !scalars
2606 : integer,intent(in) :: constraint,natom,trotter
2607 : !arrays
2608 : real(dp),intent(in) :: mass(natom),wtatcon(3,natom),xcart(3,natom,trotter)
2609 : real(dp),intent(out) :: constraint_output(2)
2610 : real(dp),intent(inout) :: forces(3,natom,trotter)
2611 : !Local variables-------------------------------
2612 : !scalars
2613 : integer :: iatom,ii,iimage
2614 : real(dp) :: af,lambda_cst,masstot,one_over_trotter,xcart_centroid,zz
2615 : !character(len=500) :: msg
2616 : !arrays
2617 : real(dp) :: force_centroid(3),lambda_com(3),mat(3,3),matinv(3,3),vec(3),weightsum(3)
2618 :
2619 : !************************************************************************
2620 :
2621 :
2622 510 : select case(constraint)
2623 :
2624 : !=== No constraint =======================================================
2625 : case(0)
2626 :
2627 242 : constraint_output(:)=zero
2628 242 : return
2629 :
2630 : !=== Linear combination of centroid coordinates ==========================
2631 : case(1)
2632 :
2633 : ! Some useful quantities
2634 182 : masstot=sum(mass)*dble(trotter)
2635 572 : weightsum(:)=sum(wtatcon,dim=2)*dble(trotter)
2636 26 : zz=zero;af=zero;force_centroid=zero
2637 156 : do iimage=1,trotter
2638 936 : do iatom=1,natom
2639 3250 : do ii=1,3
2640 2340 : force_centroid(ii)=force_centroid(ii)+forces(ii,iatom,iimage)
2641 2340 : af=af+wtatcon(ii,iatom)*forces(ii,iatom,iimage)/mass(iatom)
2642 3120 : zz=zz+wtatcon(ii,iatom)**2/mass(iatom)
2643 : end do
2644 : end do
2645 : end do
2646 104 : vec(:)=force_centroid(:)-weightsum(:)*af/zz
2647 104 : do ii=1,3
2648 312 : mat(:,ii)=-weightsum(:)*weightsum(ii)/zz
2649 104 : mat(ii,ii)=mat(ii,ii)+masstot
2650 : end do
2651 26 : call matr3inv(mat,matinv)
2652 :
2653 : !Calculation of a Lagrange multipliers:
2654 : ! lambda_cst: to apply the constraint
2655 : ! lambda_com: to maintain the position of the center of mass
2656 338 : lambda_com(:)=matmul(matinv,vec)
2657 104 : lambda_cst=(af-dot_product(weightsum,lambda_com))*dble(trotter)/zz
2658 :
2659 : !Modification of forces
2660 26 : one_over_trotter=one/dble(trotter)
2661 156 : do iimage=1,trotter
2662 936 : do iatom=1,natom
2663 3250 : do ii=1,3
2664 : forces(ii,iatom,iimage)=forces(ii,iatom,iimage) &
2665 : & -lambda_cst*wtatcon(ii,iatom)*one_over_trotter &
2666 3120 : & -lambda_com(ii)*mass(iatom)
2667 : end do
2668 : end do
2669 : end do
2670 :
2671 : !Computation of relevant outputs
2672 26 : constraint_output(:)=zero
2673 : !1-Reaction coordinate
2674 182 : do iatom=1,natom
2675 650 : do ii=1,3
2676 : xcart_centroid=zero
2677 2808 : do iimage=1,trotter
2678 2808 : xcart_centroid=xcart_centroid+xcart(ii,iatom,iimage)
2679 : end do
2680 624 : constraint_output(1)=constraint_output(1)+xcart_centroid*wtatcon(ii,iatom)
2681 : end do
2682 : end do
2683 26 : constraint_output(1)=constraint_output(1)/dble(trotter)
2684 : !2-Force on reaction coordinate
2685 294 : constraint_output(2)=-lambda_cst
2686 :
2687 : end select
2688 :
2689 : end subroutine pimd_apply_constraint
2690 : !!***
2691 :
2692 : !----------------------------------------------------------------------
2693 :
2694 : !!****f* m_pimd/pimd_mass_spring
2695 : !! NAME
2696 : !! pimd_mass_spring
2697 : !!
2698 : !! FUNCTION
2699 : !! Compute masses and spring constants for PIMD. Eventually apply a coordinate transformation.
2700 : !! Possible choices for the transformation:
2701 : !! 0: no transformation
2702 : !! 1: normal mode tranformation
2703 : !! 2: staging transformation
2704 : !!
2705 : !! INPUTS
2706 : !! inertmass(natom)=fictitious masses of atoms
2707 : !! kt=kT constant
2708 : !! natom=number of atoms
2709 : !! quantummass(natom)=true masses of atoms
2710 : !! transform=coordinate transformation:
2711 : !! 0: no tranformation
2712 : !! 1: normal mode transformation
2713 : !! 2: staging transformation
2714 : !! trotter=Trotter number
2715 : !!
2716 : !! OUTPUT
2717 : !!
2718 : !! SIDE EFFECTS
2719 : !! mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
2720 : !! spring(natom,mass_dim)=spring constants of atoms (mass_dim=1 or trotter)
2721 : !!
2722 : !! NOTES
2723 : !! Back transformation (ioption=-1) not implemented !
2724 : !!
2725 : !! SOURCE
2726 :
2727 45 : subroutine pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,transform,trotter)
2728 :
2729 : !Arguments ------------------------------------
2730 : !scalars
2731 : integer,intent(in) :: natom,transform,trotter
2732 : real(dp),intent(in) :: kt
2733 : !arrays
2734 : real(dp),intent(in) :: inertmass(natom),quantummass(natom)
2735 : real(dp),intent(out) :: mass(:,:),spring(:,:)
2736 : !Local variables-------------------------------
2737 : !scalars
2738 : integer :: iimage,kk,natom_mass,natom_spring,nmass,nspring
2739 : real(dp) :: gammasquare
2740 : character(len=500) :: msg
2741 : !arrays
2742 45 : real(dp),allocatable :: mass_temp(:,:),lambda(:)
2743 :
2744 : !************************************************************************
2745 :
2746 45 : natom_mass =size(mass ,1);nmass =size(mass ,2)
2747 45 : natom_spring=size(spring,1);nspring=size(spring,2)
2748 45 : if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
2749 0 : msg='Wrong dimensions for array mass !'
2750 0 : ABI_BUG(msg)
2751 : end if
2752 45 : if (natom/=natom_spring.or.(nspring/=1.and.nspring/=trotter)) then
2753 0 : msg='Wrong dimensions for array spring !'
2754 0 : ABI_BUG(msg)
2755 : end if
2756 :
2757 : !=== No transformation ===================================================
2758 45 : if (transform==0) then
2759 : !2nd dim of mass and spring = 1
2760 135 : mass(1:natom,1)=inertmass(1:natom)
2761 135 : spring(1:natom,1)=quantummass(1:natom)*dble(trotter)*kt*kt
2762 :
2763 : !=== Normal mode transformation ==========================================
2764 20 : else if (transform==1) then
2765 :
2766 15 : ABI_MALLOC(lambda,(trotter))
2767 5 : lambda(1)=zero; lambda(trotter)=four*dble(trotter)
2768 15 : do kk=2,trotter/2
2769 : lambda(2*kk-2)=four*dble(trotter)* &
2770 10 : & (one-cos(two_pi*dble(kk-1)/dble(trotter)))
2771 : lambda(2*kk-1)=four*dble(trotter)* &
2772 15 : & (one-cos(two_pi*dble(kk-1)/dble(trotter)))
2773 : end do
2774 :
2775 : !normal mode masses
2776 35 : do iimage=1,trotter
2777 95 : mass(:,iimage)=quantummass(:)*lambda(iimage)
2778 : end do
2779 :
2780 35 : do iimage=1,trotter
2781 95 : spring(:,iimage)=mass(:,iimage)*dble(trotter)*kt*kt
2782 : end do
2783 :
2784 : !fictitious masses
2785 15 : mass(:,1)=inertmass(:)
2786 :
2787 : !from 2 to P not changed except if adiabatic PIMD
2788 : !see Hone et al, JCP 124, 154103 (2006)
2789 : gammasquare=one !adiabaticity parameter
2790 : do iimage=2,trotter
2791 : mass(:,iimage)=mass(:,iimage)/gammasquare
2792 : end do
2793 :
2794 5 : ABI_FREE(lambda)
2795 :
2796 : !=== Staging transformation ==============================================
2797 15 : else if (transform==2) then
2798 :
2799 : !Fictitious masses
2800 45 : mass(1:natom,1)=inertmass(1:natom)
2801 15 : if (nmass>1) then
2802 90 : do iimage=2,trotter
2803 240 : mass(1:natom,iimage)=inertmass(1:natom)*dble(iimage)/dble(iimage-1)
2804 : end do
2805 : end if
2806 :
2807 : !Staging masses (mass_temp)
2808 60 : ABI_MALLOC(mass_temp,(natom,trotter))
2809 45 : mass_temp(1:natom,1)=zero
2810 15 : if (nmass>1) then
2811 90 : do iimage=2,trotter
2812 240 : mass_temp(1:natom,iimage)=quantummass(1:natom)*dble(iimage)/dble(iimage-1)
2813 : end do
2814 : end if
2815 :
2816 45 : spring(1:natom,1)=mass_temp(1:natom,1)*dble(trotter)*kt*kt
2817 15 : if (nspring>1) then
2818 90 : do iimage=2,trotter
2819 240 : spring(1:natom,iimage)=mass_temp(1:natom,iimage)*dble(trotter)*kt*kt
2820 : end do
2821 : end if
2822 15 : ABI_FREE(mass_temp)
2823 :
2824 : end if
2825 :
2826 45 : end subroutine pimd_mass_spring
2827 : !!***
2828 :
2829 0 : END MODULE m_pimd
2830 : !!***
|
