28 #include "../../base/base_uses.f90"
33 LOGICAL,
PARAMETER :: debug_this_module = .false.
36 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'al_system_dynamics'
53 SUBROUTINE al_particles(al, force_env, molecule_kind_set, molecule_set, &
54 particle_set, local_molecules, local_particles, group, vel)
56 TYPE(al_system_type),
POINTER :: al
57 TYPE(force_env_type),
POINTER :: force_env
58 TYPE(molecule_kind_type),
POINTER :: molecule_kind_set(:)
59 TYPE(molecule_type),
POINTER :: molecule_set(:)
60 TYPE(particle_type),
POINTER :: particle_set(:)
61 TYPE(distribution_1d_type),
POINTER :: local_molecules, local_particles
62 TYPE(mp_comm_type),
INTENT(IN) :: group
63 REAL(kind=
dp),
INTENT(INOUT),
OPTIONAL :: vel(:, :)
65 CHARACTER(len=*),
PARAMETER :: routinen =
'al_particles'
68 LOGICAL :: my_shell_adiabatic
69 TYPE(map_info_type),
POINTER :: map_info
71 CALL timeset(routinen, handle)
72 my_shell_adiabatic = .false.
73 map_info => al%map_info
75 IF (debug_this_module) &
76 CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel,
"INIT")
78 IF (al%tau_nh <= 0.0_dp)
THEN
79 CALL al_ou_step(0.5_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
80 particle_set, local_molecules, local_particles, vel)
81 IF (debug_this_module) &
82 CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel,
"post OU")
85 CALL al_ou_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
86 particle_set, local_molecules, local_particles, vel)
87 IF (debug_this_module) &
88 CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel,
"post 1st OU")
92 local_molecules, molecule_set, group, vel=vel)
94 CALL al_nh_quarter_step(al, map_info, set_half_step_vel_factors=.true.)
98 local_molecules, my_shell_adiabatic, vel=vel)
101 local_molecules, molecule_set, group, vel=vel)
102 IF (debug_this_module) &
103 CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel,
"post rescale_vel")
106 CALL al_nh_quarter_step(al, map_info, set_half_step_vel_factors=.false.)
109 CALL al_ou_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
110 particle_set, local_molecules, local_particles, vel)
111 IF (debug_this_module) &
112 CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel,
"post 2nd OU")
117 local_molecules, molecule_set, group, vel=vel)
119 CALL timestop(handle)
131 SUBROUTINE dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, label)
132 TYPE(molecule_kind_type),
POINTER :: molecule_kind_set(:)
133 TYPE(molecule_type),
POINTER :: molecule_set(:)
134 TYPE(distribution_1d_type),
POINTER :: local_molecules
135 TYPE(particle_type),
POINTER :: particle_set(:)
136 REAL(
dp),
OPTIONAL :: vel(:, :)
137 CHARACTER(len=*) :: label
139 INTEGER :: first_atom, ikind, imol, imol_local, &
140 ipart, last_atom, nmol_local
141 TYPE(molecule_type),
POINTER :: molecule
143 DO ikind = 1,
SIZE(molecule_kind_set)
144 nmol_local = local_molecules%n_el(ikind)
145 DO imol_local = 1, nmol_local
146 imol = local_molecules%list(ikind)%array(imol_local)
147 molecule => molecule_set(imol)
148 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
149 DO ipart = first_atom, last_atom
150 IF (
PRESENT(vel))
THEN
151 WRITE (unit=*, fmt=
'("VEL ",A20," IPART ",I6," V ",3F20.10)') trim(label), ipart, vel(:, ipart)
153 WRITE (unit=*, fmt=
'("PARTICLE_SET%VEL ",A20," IPART ",I6," V ",3F20.10)') trim(label), &
154 ipart, particle_set(ipart)%v(:)
159 END SUBROUTINE dump_vel
174 SUBROUTINE al_ou_step(step, al, force_env, map_info, molecule_kind_set, molecule_set, &
175 particle_set, local_molecules, local_particles, vel)
176 REAL(
dp),
INTENT(in) :: step
177 TYPE(al_system_type),
POINTER :: al
178 TYPE(force_env_type),
POINTER :: force_env
179 TYPE(map_info_type),
POINTER :: map_info
180 TYPE(molecule_kind_type),
POINTER :: molecule_kind_set(:)
181 TYPE(molecule_type),
POINTER :: molecule_set(:)
182 TYPE(particle_type),
POINTER :: particle_set(:)
183 TYPE(distribution_1d_type),
POINTER :: local_molecules, local_particles
184 REAL(kind=
dp),
INTENT(INOUT),
OPTIONAL :: vel(:, :)
186 INTEGER :: first_atom, i, ii, ikind, imap, imol, imol_local, ipart, iparticle_kind, &
187 iparticle_local, jj, last_atom, nmol_local, nparticle, nparticle_kind, nparticle_local
188 LOGICAL :: check, present_vel
189 REAL(kind=
dp) :: mass
190 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: w(:, :)
191 TYPE(atomic_kind_type),
POINTER :: atomic_kind
192 TYPE(molecule_type),
POINTER :: molecule
194 present_vel =
PRESENT(vel)
197 DO i = 1, al%loc_num_al
198 imap = map_info%map_index(i)
200 IF (al%tau_langevin > 0.0_dp)
THEN
201 map_info%v_scale(imap) = exp(-step*al%dt/al%tau_langevin)
202 map_info%s_kin(imap) = sqrt((al%nvt(i)%nkt/al%nvt(i)%degrees_of_freedom)*(1.0_dp - map_info%v_scale(imap)**2))
204 map_info%v_scale(imap) = 1.0_dp
205 map_info%s_kin(imap) = 0.0_dp
210 nparticle =
SIZE(particle_set)
211 nparticle_kind =
SIZE(local_particles%n_el)
212 ALLOCATE (w(3, nparticle))
214 check = (nparticle_kind <=
SIZE(local_particles%n_el) .AND. nparticle_kind <=
SIZE(local_particles%list))
216 check =
ASSOCIATED(local_particles%local_particle_set)
218 DO iparticle_kind = 1, nparticle_kind
219 nparticle_local = local_particles%n_el(iparticle_kind)
220 check = (nparticle_local <=
SIZE(local_particles%list(iparticle_kind)%array))
222 DO iparticle_local = 1, nparticle_local
223 ipart = local_particles%list(iparticle_kind)%array(iparticle_local)
224 w(1, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
225 w(2, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
226 w(3, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
233 DO ikind = 1,
SIZE(molecule_kind_set)
234 nmol_local = local_molecules%n_el(ikind)
235 DO imol_local = 1, nmol_local
236 imol = local_molecules%list(ikind)%array(imol_local)
237 molecule => molecule_set(imol)
238 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
239 DO ipart = first_atom, last_atom
241 atomic_kind => particle_set(ipart)%atomic_kind
243 IF (present_vel)
THEN
245 vel(jj, ipart) = vel(jj, ipart)*map_info%p_scale(jj, ii)%point + &
246 map_info%p_kin(jj, ii)%point/sqrt(mass)*w(jj, ipart)
250 particle_set(ipart)%v(jj) = particle_set(ipart)%v(jj)*map_info%p_scale(jj, ii)%point + &
251 map_info%p_kin(jj, ii)%point/sqrt(mass)*w(jj, ipart)
260 END SUBROUTINE al_ou_step
269 SUBROUTINE al_nh_quarter_step(al, map_info, set_half_step_vel_factors)
270 TYPE(al_system_type),
POINTER :: al
271 TYPE(map_info_type),
POINTER :: map_info
272 LOGICAL,
INTENT(in) :: set_half_step_vel_factors
275 REAL(kind=
dp) :: decay, delta_k
279 DO i = 1, al%loc_num_al
280 IF (al%nvt(i)%mass > 0.0_dp)
THEN
281 imap = map_info%map_index(i)
282 delta_k = 0.5_dp*(map_info%s_kin(imap) - al%nvt(i)%nkt)
283 al%nvt(i)%chi = al%nvt(i)%chi + 0.5_dp*al%dt*delta_k/al%nvt(i)%mass
284 IF (set_half_step_vel_factors)
THEN
285 decay = exp(-0.5_dp*al%dt*al%nvt(i)%chi)
286 map_info%v_scale(imap) = decay
289 al%nvt(i)%chi = 0.0_dp
290 IF (set_half_step_vel_factors)
THEN
291 map_info%v_scale(imap) = 1.0_dp
296 END SUBROUTINE al_nh_quarter_step
subroutine, public al_particles(al, force_env, molecule_kind_set, molecule_set, particle_set, local_molecules, local_particles, group, vel)
...
Type for the canonical sampling through velocity rescaling.
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public fix_atom_control(force_env, w)
allows for fix atom constraints
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Lumps all possible extended system variables into one type for easy access and passing.
Interface for the force calculations.
Defines the basic variable types.
integer, parameter, public dp
Interface to the message passing library MPI.
Define the molecule kind structure types and the corresponding functionality.
Define the data structure for the molecule information.
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Define the data structure for the particle information.
Utilities for thermostats.
subroutine, public vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, local_molecules, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
...
subroutine, public ke_region_particles(map_info, particle_set, molecule_kind_set, local_molecules, molecule_set, group, vel)
...