Branch data 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 : 3989 : 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 [ + + ]: 3989 : if(present(force_imgmov)) then
153 : 10 : imgmov=force_imgmov
154 : : else
155 : 3979 : imgmov=dtset%imgmov
156 : : end if
157 : :
158 : 3989 : call pimd_nullify(pimd_param)
159 : :
160 [ + + ]: 3989 : 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 : 3989 : 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 : 7968 : subroutine pimd_nullify(pimd_param)
227 : :
228 : : !Arguments ------------------------------------
229 : : !scalars
230 : : type(pimd_type),intent(inout) :: pimd_param
231 : :
232 : : !************************************************************************
233 : :
234 : 7968 : pimd_param%adpimd = 0
235 : 7968 : pimd_param%constraint = 0
236 : 7968 : pimd_param%irandom = -1
237 : 7968 : pimd_param%nnos = -1
238 : 7968 : pimd_param%ntypat = -1
239 : 7968 : pimd_param%optcell = -1
240 : 7968 : pimd_param%pitransform = -1
241 : 7968 : pimd_param%qtb_file_unit= -1
242 : 7968 : pimd_param%traj_unit = -1
243 : 7968 : pimd_param%use_qtb = 0
244 : 7968 : pimd_param%adpimd_gamma = one
245 : 7968 : pimd_param%friction = zero
246 : 7968 : pimd_param%bmass = zero
247 : 7968 : pimd_param%dtion = zero
248 : 7968 : pimd_param%frictionbar = zero
249 : 7968 : nullify(pimd_param%mdtemp)
250 : 7968 : nullify(pimd_param%pimass)
251 : 7968 : nullify(pimd_param%strtarget)
252 : 7968 : nullify(pimd_param%amu)
253 : 7968 : nullify(pimd_param%qmass)
254 : 7968 : nullify(pimd_param%typat)
255 : 7968 : nullify(pimd_param%wtatcon)
256 : :
257 : 7968 : 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 : 3979 : 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 [ + + ]: 3979 : ABI_SFREE(pimd_param%zeta_prev)
295 [ + + ]: 3979 : ABI_SFREE(pimd_param%zeta)
296 [ + + ]: 3979 : ABI_SFREE(pimd_param%zeta_next)
297 [ + + ]: 3979 : ABI_SFREE(pimd_param%dzeta)
298 : :
299 [ - + ]: 3979 : 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 [ - + ]: 3979 : 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 : 3979 : call pimd_nullify(pimd_param)
312 : :
313 : 3979 : 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 : : !!***
|
