bug fix AD_LANGEVIN with fixed atom constraints (Noam)

svn-origin-rev: 14922
This commit is contained in:
Teodoro Laino 2015-01-18 19:11:54 +00:00
parent cef3284141
commit 4c2ba344de
4 changed files with 38 additions and 24 deletions

View file

@ -733,7 +733,11 @@ CONTAINS
CPPrecondition(.NOT.ASSOCIATED(section),cp_failure_level,routineP,error,failure)
IF (.NOT. failure) THEN
CALL section_create(section,name="ad_langevin",&
description="Parameters of the adaptive-Langevin thermostat.",&
description="Parameters of the adaptive-Langevin thermostat."//&
" Known to work with NVT ensemble, but not tested with"//&
" other ensembles. Also tested with FIXED_ATOMS constraints, but"//&
" may not work with other constraints (restraints should be OK, but"//&
" haven't been well tested)",&
citations=(/Jones2011/),required=.TRUE.,error=error)
NULLIFY(keyword, subsection)

View file

@ -11,8 +11,10 @@ MODULE al_system_dynamics
USE al_system_types, ONLY: al_system_type
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE constraint_fxd, ONLY: fix_atom_control
USE distribution_1d_types, ONLY: distribution_1d_type
USE extended_system_types, ONLY: map_info_type
USE force_env_types, ONLY: force_env_type
USE kinds, ONLY: dp
USE molecule_kind_types, ONLY: molecule_kind_type
USE molecule_types_new, ONLY: get_molecule,&
@ -39,6 +41,7 @@ CONTAINS
! *****************************************************************************
!> \brief ...
!> \param al ...
!> \param force_env ...
!> \param molecule_kind_set ...
!> \param molecule_set ...
!> \param particle_set ...
@ -49,10 +52,11 @@ CONTAINS
!> \param error ...
!> \author Noam Bernstein [noamb] 02.2012
! *****************************************************************************
SUBROUTINE al_particles( al, molecule_kind_set, molecule_set, &
SUBROUTINE al_particles( al, force_env, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, group, vel, error)
TYPE(al_system_type), POINTER :: al
TYPE(force_env_type), POINTER :: force_env
TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : )
TYPE(molecule_type), POINTER :: molecule_set( : )
TYPE(particle_type), POINTER :: particle_set( : )
@ -77,13 +81,13 @@ CONTAINS
CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "INIT")
IF (al%tau_nh <= 0.0_dp) THEN
CALL al_OU_step(0.5_dp, al, map_info, molecule_kind_set, molecule_set, &
CALL al_OU_step(0.5_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, vel, error)
IF (debug_this_module) &
CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post OU")
ELSE
! quarter step of Langevin using Ornstein-Uhlenbeck
CALL al_OU_step(0.25_dp, al, map_info, molecule_kind_set, molecule_set, &
CALL al_OU_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, vel, error)
IF (debug_this_module) &
CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post 1st OU")
@ -107,7 +111,7 @@ CONTAINS
CALL al_NH_quarter_step(al, map_info, set_half_step_vel_factors=.FALSE., error=error)
! quarter step of Langevin using Ornstein-Uhlenbeck
CALL al_OU_step(0.25_dp, al, map_info, molecule_kind_set, molecule_set, &
CALL al_OU_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, vel, error)
IF (debug_this_module) &
CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post 2nd OU")
@ -163,6 +167,7 @@ CONTAINS
!> \brief ...
!> \param step ...
!> \param al ...
!> \param force_env ...
!> \param map_info ...
!> \param molecule_kind_set ...
!> \param molecule_set ...
@ -172,10 +177,11 @@ CONTAINS
!> \param vel ...
!> \param error ...
! *****************************************************************************
SUBROUTINE al_OU_step(step, al, map_info, molecule_kind_set, molecule_set, &
SUBROUTINE al_OU_step(step, al, force_env, map_info, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, vel, error)
REAL(dp), INTENT(in) :: step
TYPE(al_system_type), POINTER :: al
TYPE(force_env_type), POINTER :: force_env
TYPE(map_info_type), POINTER :: map_info
TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : )
TYPE(molecule_type), POINTER :: molecule_set( : )
@ -237,6 +243,8 @@ CONTAINS
END DO
END DO
CALL fix_atom_control(force_env,error,w)
ii = 0
DO ikind=1, SIZE(molecule_kind_set)
nmol_local = local_molecules % n_el (ikind)

View file

@ -869,10 +869,10 @@ CONTAINS
local_particles, para_env%group, shell_particle_set=shell_particle_set,&
core_particle_set=core_particle_set, error=error)
ELSE
CALL apply_thermostat_particles(thermostat_fast, molecule_kind_set, molecule_set,&
CALL apply_thermostat_particles(thermostat_fast, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, para_env%group ,error=error)
CALL apply_thermostat_particles(thermostat_slow, molecule_kind_set, molecule_set,&
CALL apply_thermostat_particles(thermostat_slow, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, para_env%group ,error=error)
END IF
@ -943,10 +943,10 @@ CONTAINS
local_particles, para_env%group, vel= tmp%vel, shell_vel= tmp%shell_vel, &
core_vel= tmp%core_vel, error=error)
ELSE
CALL apply_thermostat_particles(thermostat_slow,molecule_kind_set, molecule_set, &
CALL apply_thermostat_particles(thermostat_slow, force_env,molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, vel= tmp%vel, error=error)
CALL apply_thermostat_particles(thermostat_fast,molecule_kind_set, molecule_set, &
CALL apply_thermostat_particles(thermostat_fast, force_env,molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, vel= tmp%vel, error=error)
END IF
@ -1089,7 +1089,7 @@ CONTAINS
! Apply Thermostat over the full set of particles
IF(shell_adiabatic) THEN
CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic,&
shell_particle_set=shell_particle_set, core_particle_set=core_particle_set,&
error=error)
@ -1098,7 +1098,7 @@ CONTAINS
local_particles, para_env%group, shell_particle_set=shell_particle_set,&
core_particle_set=core_particle_set, error=error)
ELSE
CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, para_env%group ,error=error)
END IF
@ -1199,7 +1199,7 @@ CONTAINS
! Apply Thermostat over the full set of particles
IF (shell_adiabatic) THEN
CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
CALL apply_thermostat_particles(thermostat_part, force_env,molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic,&
vel= tmp%vel, shell_vel= tmp%shell_vel, core_vel= tmp%core_vel, error=error)
@ -1207,7 +1207,7 @@ CONTAINS
local_particles, para_env%group, vel= tmp%vel, shell_vel= tmp%shell_vel, &
core_vel= tmp%core_vel, error=error)
ELSE
CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
CALL apply_thermostat_particles(thermostat_part, force_env,molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, vel= tmp%vel, error=error)
END IF
@ -1377,13 +1377,13 @@ CONTAINS
! Apply Thermostat over the full set of particles
IF(simpar% ensemble /= npe_i_ensemble) THEN
IF(shell_adiabatic) THEN
CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, para_env%group,shell_adiabatic=shell_adiabatic,&
shell_particle_set=shell_particle_set, core_particle_set=core_particle_set,&
error=error)
ELSE
CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, para_env%group ,error=error)
END IF
END IF
@ -1525,11 +1525,11 @@ CONTAINS
! Apply Thermostat over the full set of particles
IF(simpar% ensemble /= npe_i_ensemble) THEN
IF (shell_adiabatic) THEN
CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
CALL apply_thermostat_particles(thermostat_part, force_env,molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic,&
vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel, error=error)
ELSE
CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
CALL apply_thermostat_particles(thermostat_part, force_env,molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, vel=tmp%vel, error=error)
END IF
END IF
@ -2396,12 +2396,12 @@ CPSourceFileRef,&
! Apply Thermostat over the full set of particles
IF(simpar% ensemble /= npe_f_ensemble) THEN
IF(shell_adiabatic) THEN
CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, para_env%group,shell_adiabatic=shell_adiabatic,&
shell_particle_set=shell_particle_set,core_particle_set=core_particle_set,&
error=error)
ELSE
CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, para_env%group ,error=error)
END IF
END IF
@ -2539,12 +2539,12 @@ CPSourceFileRef,&
! Apply Thermostat over the full set of particles
IF(simpar% ensemble /= npe_f_ensemble) THEN
IF (shell_adiabatic) THEN
CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
CALL apply_thermostat_particles(thermostat_part, force_env,molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group,shell_adiabatic=shell_adiabatic,&
vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel, error=error)
ELSE
CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
CALL apply_thermostat_particles(thermostat_part, force_env,molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, vel=tmp%vel, error=error)
END IF
END IF

View file

@ -580,6 +580,7 @@ CPSourceFileRef,&
! *****************************************************************************
!> \brief ...
!> \param thermostat ...
!> \param force_env ...
!> \param molecule_kind_set ...
!> \param molecule_set ...
!> \param particle_set ...
@ -598,11 +599,12 @@ CPSourceFileRef,&
!> 10.2007 created [tlaino]
!> \author Teodoro Laino
! *****************************************************************************
SUBROUTINE apply_thermostat_particles(thermostat, molecule_kind_set, molecule_set,&
SUBROUTINE apply_thermostat_particles(thermostat, force_env, molecule_kind_set, molecule_set,&
particle_set, local_molecules, local_particles, group, shell_adiabatic, shell_particle_set,&
core_particle_set, vel, shell_vel, core_vel, error)
TYPE(thermostat_type), POINTER :: thermostat
TYPE(force_env_type), POINTER :: force_env
TYPE(molecule_kind_type), POINTER :: molecule_kind_set( : )
TYPE(molecule_type), POINTER :: molecule_set( : )
TYPE(particle_type), POINTER :: particle_set( : )
@ -638,7 +640,7 @@ CPSourceFileRef,&
ELSE IF (thermostat%type_of_thermostat==do_thermo_al) THEN
! Apply AD_LANGEVIN Thermostat
CPPostcondition(ASSOCIATED(thermostat%al),cp_failure_level,routineP,error,failure)
CALL al_particles(thermostat%al,molecule_kind_set,molecule_set,&
CALL al_particles(thermostat%al,force_env,molecule_kind_set,molecule_set,&
particle_set,local_molecules, local_particles, group, vel, error)
ELSE IF (thermostat%type_of_thermostat==do_thermo_gle) THEN
! Apply GLE Thermostat