comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

VIRTUAL SITE constraint. This allows for TIP4P water. Added regtest (patch

Browse source 
created by Marcel Baer and fixed/cleaned by Teo). Reset few regtests for numerics.


svn-origin-rev: 8034
This commit is contained in:
Teodoro Laino 2008-12-09 16:18:01 +00:00
parent 93e8156b48
commit ee457522e3
19 changed files with 1307 additions and 288 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/OBJECTDEFS
View file

@ -46,6 +46,7 @@ OBJECTS_GENERIC =\
constraint_clv.o\
constraint_fxd.o\
constraint_util.o\
constraint_vsite.o\
core_ppnl.o\
cp2k_debug.o\
cp2k_info.o\

46
src/constraint.F
View file

@ -37,6 +37,8 @@ MODULE constraint
get_roll_matrix,&
restore_temporary_set,&
update_temporary_set
USE constraint_vsite, ONLY: shake_vsite_ext,&
shake_vsite_int
USE cp_para_types, ONLY: cp_para_env_type
USE distribution_1d_types, ONLY: distribution_1d_type
USE f77_blas
@ -100,7 +102,7 @@ CONTAINS
routineP = moduleN//':'//routineN
INTEGER :: handle, i, ikind, imol, ishake_ext, ishake_int, k, n3x3con, &
n4x6con, nconstraint, nkind, nmol_per_kind
n4x6con, nconstraint, nkind, nmol_per_kind, nvsitecon
LOGICAL :: do_ext_constraint
REAL(KIND=dp) :: int_max_sigma, mass, max_sigma
REAL(KIND=dp), DIMENSION(SIZE(pos, 2)) :: imass
@ -114,7 +116,7 @@ CONTAINS
DO k = 1, SIZE(pos,2)
atomic_kind=>particle_set (k) % atomic_kind
CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass)
imass(k) = 1.0_dp/mass
imass(k) = 1.0_dp/MAX(mass,EPSILON(0.0_dp))
END DO
do_ext_constraint = (gci%ntot/=0)
ishake_ext = 0
@ -131,7 +133,7 @@ CONTAINS
molecule_kind => molecule % molecule_kind
CALL get_molecule_kind ( molecule_kind, ncolv = ncolv, &
ng3x3 = n3x3con, ng4x6=n4x6con,&
nconstraint=nconstraint )
nconstraint=nconstraint,nvsite=nvsitecon)
IF (nconstraint==0) CYCLE
ishake_int = 0
int_max_sigma = -1.0E+10_dp
@ -153,6 +155,10 @@ CONTAINS
END DO Shake_Intra_Loop
max_sigma = MAX(max_sigma, int_max_sigma)
CALL shake_int_info(log_unit,i,ishake_int,max_sigma,error)
! Virtual Site
IF ( nvsitecon /= 0 ) &
CALL shake_vsite_int( molecule, particle_set, pos, vel, dt, ishake_int,&
int_max_sigma, error)
END DO
END DO MOL
! Intermolecular constraints
@ -170,12 +176,17 @@ CONTAINS
IF ( gci%ncolv%ntot /= 0 ) &
CALL shake_colv_ext(gci, particle_set, pos, vel, dt, ishake_ext,&
cell, imass, max_sigma, error)
! Virtual Site
IF ( gci%nvsite /= 0 ) &
CALL shake_vsite_ext( gci, particle_set, pos, vel, dt, ishake_ext,&
max_sigma, error)
CALL restore_temporary_set( particle_set, local_particles, pos=pos, vel=vel )
END IF
CALL shake_ext_info(log_unit,ishake_ext,max_sigma,error)
END DO Shake_Inter_Loop
CALL dump_lagrange_mult(dump_lm, lagrange_mult,local_molecules,molecule_set,gci,&
molecule_kind_set,group,"S",error)
CALL timestop(handle)
END SUBROUTINE shake_control
@ -219,7 +230,7 @@ CONTAINS
DO k = 1, SIZE(vel,2)
atomic_kind=>particle_set (k) % atomic_kind
CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass)
imass(k) = 1.0_dp/mass
imass(k) = 1.0_dp/MAX(mass,EPSILON(0.0_dp))
END DO
do_ext_constraint = (gci%ntot/=0)
irattle_ext = 0
@ -311,7 +322,7 @@ CONTAINS
INTEGER :: handle, i, ikind, imol, ishake_ext, ishake_int, k, &
lagrange_mult, log_unit, n3x3con, n4x6con, nconstraint, nkind, &
nmol_per_kind
nmol_per_kind, nvsitecon
LOGICAL :: do_ext_constraint, dump_lm
REAL(KIND=dp) :: int_max_sigma, mass, &
max_sigma, shake_tol
@ -337,7 +348,7 @@ CONTAINS
DO k = 1, SIZE(pos,2)
atomic_kind=>particle_set (k) % atomic_kind
CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass)
imass(k) = 1.0_dp/mass
imass(k) = 1.0_dp/MAX(mass,EPSILON(0.0_dp))
END DO
do_ext_constraint = (gci%ntot/=0)
ishake_ext = 0
@ -354,7 +365,7 @@ CONTAINS
molecule_kind => molecule % molecule_kind
CALL get_molecule_kind ( molecule_kind, ncolv = ncolv, &
ng3x3 = n3x3con, ng4x6=n4x6con,&
nconstraint=nconstraint )
nconstraint=nconstraint,nvsite=nvsitecon)
IF (nconstraint==0) CYCLE
ishake_int = 0
int_max_sigma = -1.0E+10_dp
@ -376,6 +387,12 @@ CONTAINS
END DO Shake_Roll_Intra_Loop
max_sigma = MAX(max_sigma, int_max_sigma)
CALL shake_int_info(log_unit,i,ishake_int,max_sigma,error)
! Virtual Site
IF ( nvsitecon /= 0 ) THEN
CALL cp_unimplemented_error(fromWhere=routineP, &
message="Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!", &
error=error, error_level=cp_failure_level)
END IF
END DO
END DO MOL
! Intermolecular constraints
@ -393,6 +410,11 @@ CONTAINS
IF ( gci%ncolv%ntot /= 0 ) &
CALL shake_roll_colv_ext( gci, particle_set, pos, vel, r_shake,&
v_shake, dt, ishake_ext, cell, imass, max_sigma, error )
! Virtual Site
IF ( gci%nvsite /= 0 ) &
CALL cp_unimplemented_error(fromWhere=routineP, &
message="Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!", &
error=error, error_level=cp_failure_level)
CALL restore_temporary_set ( particle_set, local_particles, pos = pos, vel = vel )
END IF
CALL shake_ext_info(log_unit,ishake_ext,max_sigma,error)
@ -410,7 +432,7 @@ CONTAINS
! *****************************************************************************
SUBROUTINE rattle_roll_control( gci, local_molecules, molecule_set, &
molecule_kind_set, particle_set, vel, dt, simpar, virial, vector,&
veps, mass, roll_tol, iroll, para_env, u, cell, local_particles, error )
veps, roll_tol, iroll, para_env, u, cell, local_particles, error )
TYPE(global_constraint_type), POINTER :: gci
TYPE(distribution_1d_type), POINTER :: local_molecules
@ -424,8 +446,6 @@ CONTAINS
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vector
REAL(KIND=dp), DIMENSION(:, :), &
INTENT(INOUT) :: veps
REAL(KIND=dp), DIMENSION(:, :), &
INTENT(IN) :: mass
REAL(KIND=dp), INTENT(OUT) :: roll_tol
INTEGER, INTENT(INOUT) :: iroll
TYPE(cp_para_env_type), INTENT(IN) :: para_env
@ -442,7 +462,7 @@ CONTAINS
lagrange_mult, log_unit, n3x3con, n4x6con, nconstraint, nkind, &
nmol_per_kind
LOGICAL :: do_ext_constraint, dump_lm
REAL(KIND=dp) :: amass, int_max_sigma, &
REAL(KIND=dp) :: int_max_sigma, mass, &
max_sigma, rattle_tol
REAL(KIND=dp), DIMENSION(3, 3) :: r_rattle
REAL(KIND=dp), DIMENSION(SIZE(vel, 2)) :: imass
@ -466,8 +486,8 @@ CONTAINS
END IF
DO k = 1, SIZE(vel,2)
atomic_kind=>particle_set (k) % atomic_kind
CALL get_atomic_kind(atomic_kind=atomic_kind,mass=amass)
imass(k) = 1.0_dp/amass
CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass)
imass(k) = 1.0_dp/MAX(mass,EPSILON(0.0_dp))
END DO
do_ext_constraint = (gci%ntot/=0)
irattle_ext = 0

293
src/constraint_vsite.F Normal file
View file

@ -0,0 +1,293 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2008 CP2K developers group !
!-----------------------------------------------------------------------------!
! *****************************************************************************
!> \brief Routines to handle the virtual site constraint/restraint
!> \par History
!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
!> (patch by Marcel Baer)
! *****************************************************************************
MODULE constraint_vsite
USE cp_subsys_types, ONLY: cp_subsys_get,&
cp_subsys_type
USE distribution_1d_types, ONLY: distribution_1d_type
USE f77_blas
USE force_env_types, ONLY: force_env_get,&
force_env_type
USE kinds, ONLY: dp
USE mol_kind_new_list_types, ONLY: mol_kind_new_list_type
USE mol_new_list_types, ONLY: mol_new_list_type
USE molecule_kind_types, ONLY: fixd_constraint_type,&
get_molecule_kind,&
molecule_kind_type,&
vsite_constraint_type
USE molecule_types_new, ONLY: get_molecule,&
global_constraint_type,&
molecule_type
USE particle_list_types, ONLY: particle_list_type
USE particle_types, ONLY: particle_type
#include "cp_common_uses.h"
IMPLICIT NONE
PRIVATE
PUBLIC :: shake_vsite_int,&
force_vsite_int,&
force_vsite_ext,&
shake_vsite_ext,&
vsite_force_control
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_vsite'
CONTAINS
! *****************************************************************************
!> \brief control force distribution for virtual sites
!> \par History
!> - none
!> \author Marcel Baer
!> \date 12.2008
! *****************************************************************************
SUBROUTINE vsite_force_control( force_env, error,w)
TYPE(force_env_type), POINTER :: force_env
TYPE(cp_error_type), INTENT(inout) :: error
REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: w
INTEGER :: i, ikind, imol, nconstraint, &
nkind, nmol_per_kind, &
nvsitecon
LOGICAL :: do_ext_constraint
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(distribution_1d_type), POINTER :: local_molecules, &
local_particles
TYPE(global_constraint_type), POINTER :: gci
TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds
TYPE(mol_new_list_type), POINTER :: molecules
TYPE(molecule_kind_type), DIMENSION(:), &
POINTER :: molecule_kind_set
TYPE(molecule_kind_type), POINTER :: molecule_kind
TYPE(molecule_type), DIMENSION(:), &
POINTER :: molecule_set
TYPE(molecule_type), POINTER :: molecule
TYPE(particle_list_type), POINTER :: particles
TYPE(particle_type), DIMENSION(:), &
POINTER :: particle_set
NULLIFY(gci, subsys, local_molecules, local_particles,&
molecule_kinds)
CALL force_env_get(force_env=force_env, subsys=subsys,error=error)
CALL cp_subsys_get(subsys=subsys,local_particles=local_particles,&
particles=particles,local_molecules_new=local_molecules,&
molecule_kinds_new=molecule_kinds,gci=gci,molecules_new=molecules,&
error=error)
molecule_kind_set => molecule_kinds%els
molecule_set => molecules%els
particle_set => particles%els
nkind = SIZE ( molecule_kind_set )
! Intermolecular Virtual Site Constraints
do_ext_constraint = .FALSE.
IF (ASSOCIATED(gci)) THEN
do_ext_constraint = (gci%ntot/=0)
END IF
! Intramolecular Virtual Site Constraints
MOL: DO ikind = 1, nkind
nmol_per_kind = local_molecules % n_el ( ikind )
DO imol = 1, nmol_per_kind
i = local_molecules % list ( ikind ) % array ( imol )
molecule => molecule_set ( i )
molecule_kind => molecule % molecule_kind
CALL get_molecule_kind ( molecule_kind,nconstraint=nconstraint,nvsite=nvsitecon)
IF (nconstraint==0) CYCLE
IF (nvsitecon/=0) &
CALL force_vsite_int( molecule, particle_set)
END DO
END DO MOL
! Intermolecular Virtual Site Constraints
IF (do_ext_constraint) THEN
IF ( gci%nvsite /= 0 ) &
CALL force_vsite_ext( gci, particle_set)
END IF
END SUBROUTINE vsite_force_control
! *****************************************************************************
!> \brief Intramolecular virtual site
!> \par History
!> 12.2008 Marcel Baer
! *****************************************************************************
SUBROUTINE shake_vsite_int( molecule, particle_set, pos, vel, dt, ishake,&
max_sigma, error)
TYPE(molecule_type), POINTER :: molecule
TYPE(particle_type), POINTER :: particle_set( : )
REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : )
REAL(kind=dp), INTENT(in) :: dt
INTEGER, INTENT(IN) :: ishake
REAL(KIND=dp), INTENT(INOUT) :: max_sigma
TYPE(cp_error_type), INTENT(inout) :: error
INTEGER :: first_atom, nvsite
TYPE(fixd_constraint_type), &
DIMENSION(:), POINTER :: fixd_list
TYPE(molecule_kind_type), POINTER :: molecule_kind
TYPE(vsite_constraint_type), POINTER :: vsite_list( : )
NULLIFY(fixd_list)
molecule_kind => molecule % molecule_kind
CALL get_molecule_kind ( molecule_kind, nvsite = nvsite,&
fixd_list=fixd_list, vsite_list = vsite_list )
CALL get_molecule ( molecule, first_atom = first_atom)
! Real Shake
CALL shake_vsite_low( fixd_list, vsite_list, nvsite, first_atom, &
particle_set, pos, vel, dt, ishake, max_sigma, error)
END SUBROUTINE shake_vsite_int
! *****************************************************************************
!> \brief Intramolecular virtual site
!> \par History
!> 12.2008 Marcel Baer
! *****************************************************************************
SUBROUTINE shake_vsite_ext( gci, particle_set, pos, vel, dt, ishake,&
max_sigma, error)
TYPE(global_constraint_type), POINTER :: gci
TYPE(particle_type), POINTER :: particle_set( : )
REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : )
REAL(kind=dp), INTENT(in) :: dt
INTEGER, INTENT(IN) :: ishake
REAL(KIND=dp), INTENT(INOUT) :: max_sigma
TYPE(cp_error_type), INTENT(inout) :: error
INTEGER :: first_atom, nvsite
TYPE(fixd_constraint_type), &
DIMENSION(:), POINTER :: fixd_list
TYPE(vsite_constraint_type), POINTER :: vsite_list( : )
first_atom = 1
nvsite = gci%nvsite
fixd_list => gci%fixd_list
vsite_list => gci%vsite_list
! Real Shake
CALL shake_vsite_low( fixd_list, vsite_list, nvsite, first_atom, &
particle_set, pos, vel, dt, ishake, max_sigma, error)
END SUBROUTINE shake_vsite_ext
! *****************************************************************************
!> \par History
!> 12.2008 Marcel Bear
! *****************************************************************************
SUBROUTINE shake_vsite_low(fixd_list, vsite_list, nvsite, first_atom, &
particle_set, pos, vel, dt, ishake, max_sigma, error)
TYPE(fixd_constraint_type), &
DIMENSION(:), POINTER :: fixd_list
TYPE(vsite_constraint_type), POINTER :: vsite_list( : )
INTEGER, INTENT(IN) :: nvsite, first_atom
TYPE(particle_type), POINTER :: particle_set( : )
REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : )
REAL(kind=dp), INTENT(in) :: dt
INTEGER, INTENT(IN) :: ishake
REAL(KIND=dp), INTENT(INOUT) :: max_sigma
TYPE(cp_error_type), INTENT(inout) :: error
INTEGER :: iconst, index_a, index_b, &
index_c, index_d
REAL(KIND=dp), DIMENSION(3) :: r1, r2
DO iconst = 1, nvsite
IF (vsite_list(iconst)%restraint%active) CYCLE
index_a = vsite_list(iconst)%a + first_atom -1
index_b = vsite_list(iconst)%b + first_atom -1
index_c = vsite_list(iconst)%c + first_atom -1
index_d = vsite_list(iconst)%d + first_atom -1
r1(:) = pos(:,index_b) - pos(:,index_c )
r2(:) = pos(:,index_d) - pos(:,index_c )
pos(:,index_a)=pos(:,index_c)+vsite_list(iconst)%wbc*r1(:)+&
vsite_list(iconst)%wdc*r2(:)
END DO
END SUBROUTINE shake_vsite_low
! *****************************************************************************
!> \brief Intramolecular virtual site
!> \par History
!> 12.2008 Marcel Bear
! *****************************************************************************
SUBROUTINE force_vsite_int( molecule, particle_set)
TYPE(molecule_type), POINTER :: molecule
TYPE(particle_type), POINTER :: particle_set( : )
INTEGER :: first_atom, iconst, index_a, &
index_b, index_c, index_d, &
nvsite
REAL(KIND=dp) :: wb, wc, wd
TYPE(molecule_kind_type), POINTER :: molecule_kind
TYPE(vsite_constraint_type), POINTER :: vsite_list( : )
molecule_kind => molecule % molecule_kind
CALL get_molecule_kind ( molecule_kind, nvsite = nvsite, vsite_list = vsite_list )
CALL get_molecule ( molecule, first_atom = first_atom)
DO iconst = 1, nvsite
IF (vsite_list(iconst)%restraint%active) CYCLE
index_a = vsite_list(iconst)%a + first_atom -1
index_b = vsite_list(iconst)%b + first_atom -1
index_c = vsite_list(iconst)%c + first_atom -1
index_d = vsite_list(iconst)%d + first_atom -1
wb=vsite_list(iconst)%wbc
wd=vsite_list(iconst)%wdc
wc=1.0_dp-vsite_list(iconst)%wbc-vsite_list(iconst)%wdc
particle_set(index_b)%f(:)=particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
particle_set(index_c)%f(:)=particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
particle_set(index_d)%f(:)=particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
particle_set(index_a)%f(:)=0.0_dp
END DO
END SUBROUTINE force_vsite_int
! *****************************************************************************
!> \brief Intramolecular virtual site
!> \par History
!> 12.2008 Marcel Bear
! *****************************************************************************
SUBROUTINE force_vsite_ext( gci, particle_set)
TYPE(global_constraint_type), POINTER :: gci
TYPE(particle_type), POINTER :: particle_set( : )
INTEGER :: first_atom, iconst, index_a, &
index_b, index_c, index_d, &
nvsite
REAL(KIND=dp) :: wb, wc, wd
TYPE(vsite_constraint_type), POINTER :: vsite_list( : )
first_atom = 1
nvsite = gci%nvsite
vsite_list => gci%vsite_list
! Real Shake
DO iconst = 1, nvsite
IF (vsite_list(iconst)%restraint%active) CYCLE
index_a = vsite_list(iconst)%a + first_atom -1
index_b = vsite_list(iconst)%b + first_atom -1
index_c = vsite_list(iconst)%c + first_atom -1
index_d = vsite_list(iconst)%d + first_atom -1
wb=vsite_list(iconst)%wbc
wd=vsite_list(iconst)%wdc
wc=1.0_dp-vsite_list(iconst)%wbc-vsite_list(iconst)%wdc
particle_set(index_b)%f(:)=particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
particle_set(index_c)%f(:)=particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
particle_set(index_d)%f(:)=particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
particle_set(index_a)%f(:)=0.0_dp
END DO
END SUBROUTINE force_vsite_ext
END MODULE constraint_vsite

4
src/force_env_methods.F
View file

@ -20,6 +20,7 @@ MODULE force_env_methods
shake_control
USE constraint_fxd, ONLY: fix_atom_control
USE constraint_util, ONLY: getold
USE constraint_vsite, ONLY: vsite_force_control
USE cp_output_handling, ONLY: cp_p_file,&
cp_print_key_finished_output,&
cp_print_key_should_output,&
@ -299,6 +300,8 @@ CONTAINS
CALL fix_atom_control(force_env, error=error)
! All Restraints
CALL restraint_control(force_env, error=error)
! Virtual Sites
CALL vsite_force_control(force_env,error)
! External Potential
CALL add_external_potential(force_env, error=error)
! Rescale forces if requested
@ -1488,7 +1491,6 @@ CPSourceFileRef,&
TYPE(particle_list_type), POINTER :: particles
failure=.FALSE.
CALL timeset(routineN,handle)
CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure)
CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure)

91
src/input_cp2k_constraints.F
View file

@ -6,6 +6,8 @@
! *****************************************************************************
!> \par History
!> 10.2005 split input_cp2k into smaller modules [fawzi]
!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
!> (patch by Marcel Baer)
!> \author teo & fawzi
! *****************************************************************************
MODULE input_cp2k_constraints
@ -106,6 +108,11 @@ CONTAINS
CALL restraint_info_section(subsection, error=error)
CALL section_add_subsection(section, subsection, error=error)
CALL section_release(subsection,error=error)
CALL create_vsite_section(subsection, error=error)
CALL restraint_info_section(subsection, error=error)
CALL section_add_subsection(section, subsection, error=error)
CALL section_release(subsection,error=error)
CALL create_collective_section(subsection, error=error)
CALL restraint_info_section(subsection, error=error)
@ -708,4 +715,88 @@ CONTAINS
END SUBROUTINE create_g4x6_section
! *****************************************************************************
!> \brief Create the constraint section specialized on vsite constraints
!> \param section the section to create
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
!> \author marcel baer
! *****************************************************************************
SUBROUTINE create_vsite_section(section,error)
TYPE(section_type), POINTER :: section
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'create_vsite_section', &
routineP = moduleN//':'//routineN
LOGICAL :: failure
TYPE(keyword_type), POINTER :: keyword
failure=.FALSE.
CPPrecondition(.NOT.ASSOCIATED(section),cp_failure_level,routineP,error,failure)
IF (.NOT. failure) THEN
CALL section_create(section,name="virtual_site",&
description="This section is used to set a virtual interaction-site constraint.",&
n_keywords=3, n_subsections=0, repeats=.TRUE., required=.TRUE.,&
error=error)
NULLIFY(keyword)
! Integer
CALL keyword_create(keyword, name="MOLECULE",&
variants=(/"MOL"/),&
description="Specifies the molecule number on which constraint will be applied."//&
" MOLECULE and MOLNAME keyword exclude themself mutually.",&
usage="MOL {integer}", required=.TRUE., n_var=1, type_of_var=integer_t, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="MOLNAME",&
variants=(/"SEGNAME"/),&
description="Specifies the name of the molecule on which the constraint will be applied.",&
usage="MOLNAME {character}", required=.TRUE., n_var=1, type_of_var=char_t, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="INTERMOLECULAR",&
description="Specify if the constraint/restraint is intermolecular.",&
usage="INTERMOLECULAR <LOGICAL>",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="ATOMS",&
description="Atoms' index on which apply the constraint (v i j k), first is virtual site",&
usage="ATOMS 1 2 3 4",&
n_var=4,type_of_var=integer_t,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
! Real
CALL keyword_create(keyword, name="PARAMETERS",&
description="The constrained paramters' values to construct virtual site."//&
"r_v=a*r_ij+b*r_kj",&
usage="PARAMETERS {real} {real}",&
type_of_var=real_t, n_var=2,unit_str="internal_cp2k",error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
! Logical
CALL keyword_create(keyword, name="EXCLUDE_QM",&
description="Does not apply the constraint to the QM region within a QM/MM calculation",&
usage="EXCLUDE_QM <LOGICAL>",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="EXCLUDE_MM",&
description="Does not apply the constraint to the MM region within a QM/MM calculation",&
usage="EXCLUDE_MM <LOGICAL>",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
END IF
END SUBROUTINE create_vsite_section
END MODULE input_cp2k_constraints

464
src/integrator_utils.F
View file

@ -7,6 +7,8 @@
!> \brief Provides integrator utility routines for the integrators
!> \par History
!> Teodoro Laino [tlaino] 02.2008 - Splitting from integrator and creation
!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
!> (patch by Marcel Baer)
! *****************************************************************************
MODULE integrator_utils
@ -375,21 +377,23 @@ CONTAINS
DO iparticle_kind=1,nparticle_kind
atomic_kind => atomic_kind_set(iparticle_kind)
CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass)
nparticle_local = local_particles%n_el(iparticle_kind)
IF(my_tmpv) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
K=K+ 0.5_dp*mass *DOT_PRODUCT(tmp%vel(:,iparticle),tmp%vel(:,iparticle))
a=a+ DOT_PRODUCT(tmp%vel(:,iparticle),particle_set(iparticle)%f(:))
b=b+(1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:),particle_set(iparticle)%f(:))
END DO
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
K=K+ 0.5_dp*mass *DOT_PRODUCT(particle_set(iparticle)%v(:),particle_set(iparticle)%v(:))
a=a+ DOT_PRODUCT(particle_set(iparticle)%v(:),particle_set(iparticle)%f(:))
b=b+(1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:),particle_set(iparticle)%f(:))
END DO
IF (mass/=0.0_dp) THEN
nparticle_local = local_particles%n_el(iparticle_kind)
IF(my_tmpv) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
K=K+ 0.5_dp*mass *DOT_PRODUCT(tmp%vel(:,iparticle),tmp%vel(:,iparticle))
a=a+ DOT_PRODUCT(tmp%vel(:,iparticle),particle_set(iparticle)%f(:))
b=b+(1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:),particle_set(iparticle)%f(:))
END DO
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
K=K+ 0.5_dp*mass *DOT_PRODUCT(particle_set(iparticle)%v(:),particle_set(iparticle)%v(:))
a=a+ DOT_PRODUCT(particle_set(iparticle)%v(:),particle_set(iparticle)%f(:))
b=b+(1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:),particle_set(iparticle)%f(:))
END DO
END IF
END IF
END DO
CALL mp_sum(K,para_env % group)
@ -746,8 +750,7 @@ CONTAINS
! resetting
CALL set ( old, atomic_kind_set, particle_set, vel, local_particles, cell , npt, 'B' )
CALL rattle_roll_control( gci, local_molecules, molecule_set, molecule_kind_set, &
particle_set, vel, dt, simpar, virial, vector_v, &
npt_loc % v, npt_loc % mass, roll_tol, iroll, &
particle_set, vel, dt, simpar, virial, vector_v, npt_loc%v, roll_tol, iroll,&
para_env, u, cell, local_particles, error=error )
END SUBROUTINE rattle_roll_setup
@ -1011,140 +1014,154 @@ CONTAINS
DO iparticle_kind=1,nparticle_kind
atomic_kind => atomic_kind_set(iparticle_kind)
CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass,shell_active=is_shell)
dm = 0.5_dp*dt/mass
IF(is_shell .AND. shell_adiabatic) THEN
CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
dms = 0.5_dp*dt/shell%mass_shell
dmc = 0.5_dp*dt/shell%mass_core
fac_masss = shell%mass_shell/mass
fac_massc = shell%mass_core/mass
nparticle_local = local_particles%n_el(iparticle_kind)
IF(PRESENT(u)) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
! Transform positions and velocities and forces of the shell
CALL transform_first(shell_particle_set,tmp%shell_pos,tmp%shell_vel,&
shell_index,u,dms,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r)
IF (mass/=0.0_dp) THEN
dm = 0.5_dp*dt/mass
IF(is_shell .AND. shell_adiabatic) THEN
CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
dms = 0.5_dp*dt/shell%mass_shell
dmc = 0.5_dp*dt/shell%mass_core
fac_masss = shell%mass_shell/mass
fac_massc = shell%mass_core/mass
nparticle_local = local_particles%n_el(iparticle_kind)
IF(PRESENT(u)) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
! Transform positions and velocities and forces of the shell
CALL transform_first(shell_particle_set,tmp%shell_pos,tmp%shell_vel,&
shell_index,u,dms,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r)
! Transform positions and velocities and forces of the core
CALL transform_first(core_particle_set,tmp%core_pos,tmp%core_vel,&
shell_index,u,dmc,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r)
! Transform positions and velocities and forces of the core
CALL transform_first(core_particle_set,tmp%core_pos,tmp%core_vel,&
shell_index,u,dmc,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r)
! Derive velocities and positions of the COM
tmp%vel (:,iparticle) = fac_masss*tmp%shell_vel (:,shell_index)+&
fac_massc*tmp%core_vel (:,shell_index)
tmp%pos (:,iparticle) = fac_masss*tmp%shell_pos (:,shell_index)+&
fac_massc*tmp%core_pos (:,shell_index)
! Derive velocities and positions of the COM
tmp%vel (:,iparticle) = fac_masss*tmp%shell_vel (:,shell_index)+&
fac_massc*tmp%core_vel (:,shell_index)
tmp%pos (:,iparticle) = fac_masss*tmp%shell_pos (:,shell_index)+&
fac_massc*tmp%core_pos (:,shell_index)
tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel(1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle)))
tmp%max_vel_sc = MAX(tmp%max_vel_sc,ABS(tmp%shell_vel(1,shell_index)-tmp%core_vel(1,shell_index)),&
ABS(tmp%shell_vel(2,shell_index)-tmp%core_vel(2,shell_index)),&
ABS(tmp%shell_vel(3,shell_index)-tmp%core_vel(3,shell_index)))
tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),&
ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),&
ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel,ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),&
ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),&
ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle)))
tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel(1,iparticle)),&
ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle)))
tmp%max_vel_sc = MAX(tmp%max_vel_sc,&
ABS(tmp%shell_vel(1,shell_index)-tmp%core_vel(1,shell_index)),&
ABS(tmp%shell_vel(2,shell_index)-tmp%core_vel(2,shell_index)),&
ABS(tmp%shell_vel(3,shell_index)-tmp%core_vel(3,shell_index)))
tmp%max_dr = MAX(tmp%max_dr,&
ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),&
ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),&
ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel,&
ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),&
ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),&
ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle)))
dsc(:) = tmp%shell_pos(:,shell_index)-tmp%core_pos(:,shell_index)- &
shell_particle_set(shell_index)%r(:)+core_particle_set(shell_index)%r(:)
tmp%max_dsc = MAX(tmp%max_dsc,ABS(dsc(1)),ABS(dsc(2)),ABS(dsc(3)))
dvsc(:) = tmp%shell_vel(:,shell_index)-tmp%core_vel(:,shell_index)- &
shell_particle_set(shell_index)%v(:)+core_particle_set(shell_index)%v(:)
tmp%max_dvel_sc = MAX(tmp%max_dvel_sc,ABS(dvsc(1)),ABS(dvsc(2)),ABS(dvsc(3)))
END DO ! iparticle_local
dsc(:) = tmp%shell_pos(:,shell_index)-tmp%core_pos(:,shell_index)- &
shell_particle_set(shell_index)%r(:)+core_particle_set(shell_index)%r(:)
tmp%max_dsc = MAX(tmp%max_dsc,ABS(dsc(1)),ABS(dsc(2)),ABS(dsc(3)))
dvsc(:) = tmp%shell_vel(:,shell_index)-tmp%core_vel(:,shell_index)- &
shell_particle_set(shell_index)%v(:)+core_particle_set(shell_index)%v(:)
tmp%max_dvel_sc = MAX(tmp%max_dvel_sc,ABS(dvsc(1)),ABS(dvsc(2)),ABS(dvsc(3)))
END DO ! iparticle_local
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
tmp%shell_vel (:,shell_index) = &
shell_particle_set(shell_index)%v(:) * tmp%scale_v(:) * tmp%scale_v(:) +&
tmp%scale_v(:) * tmp%poly_v(:) * dms * shell_particle_set(shell_index)%f(:)
tmp%shell_pos (:,shell_index) = &
shell_particle_set(shell_index)%r(:) * tmp%scale_r(:) * tmp%scale_r(:) +&
tmp%scale_r(:) * tmp%poly_r(:) *dt * tmp%shell_vel ( :, shell_index )
tmp%core_vel (:,shell_index) = &
core_particle_set(shell_index)%v(:) * tmp%scale_v(:) * tmp%scale_v(:) +&
tmp%scale_v(:) * tmp%poly_v(:) *dmc * core_particle_set(shell_index)%f(:)
tmp%core_pos (:,shell_index) = &
core_particle_set(shell_index)%r(:) * tmp%scale_r(:) * tmp%scale_r(:) +&
tmp%scale_r(:) * tmp%poly_r(:) *dt * tmp%core_vel ( :, shell_index )
tmp%vel (:,iparticle) = fac_masss*tmp%shell_vel (:,shell_index)+&
fac_massc*tmp%core_vel (:,shell_index)
tmp%pos (:,iparticle) = fac_masss*tmp%shell_pos (:,shell_index)+&
fac_massc*tmp%core_pos (:,shell_index)
tmp%max_vel = MAX(tmp%max_vel,&
ABS(tmp%vel(1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle)))
tmp%max_vel_sc = MAX(tmp%max_vel_sc,&
ABS(tmp%shell_vel(1,shell_index)-tmp%core_vel(1,shell_index)),&
ABS(tmp%shell_vel(2,shell_index)-tmp%core_vel(2,shell_index)),&
ABS(tmp%shell_vel(3,shell_index)-tmp%core_vel(3,shell_index)))
tmp%max_dr = MAX(tmp%max_dr,&
ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),&
ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),&
ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel,&
ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),&
ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),&
ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle)))
dsc(:) = tmp%shell_pos (:,shell_index)-tmp%core_pos (:,shell_index)-&
shell_particle_set(shell_index)%r(:)+core_particle_set(shell_index)%r(:)
tmp%max_dsc = MAX(tmp%max_dsc,ABS(dsc(1)),ABS(dsc(2)),ABS(dsc(3)))
dvsc(:) = tmp%shell_vel(:,shell_index)-tmp%core_vel(:,shell_index)- &
shell_particle_set(shell_index)%v(:)+core_particle_set(shell_index)%v(:)
tmp%max_dvel_sc = MAX(tmp%max_dvel_sc,ABS(dvsc(1)),ABS(dvsc(2)),ABS(dvsc(3)))
END DO ! iparticle_local
END IF
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
tmp%shell_vel (:,shell_index) = &
shell_particle_set(shell_index)%v(:) * tmp%scale_v(:) * tmp%scale_v(:) +&
tmp%scale_v(:) * tmp%poly_v(:) * dms * shell_particle_set(shell_index)%f(:)
tmp%shell_pos (:,shell_index) = &
shell_particle_set(shell_index)%r(:) * tmp%scale_r(:) * tmp%scale_r(:) +&
tmp%scale_r(:) * tmp%poly_r(:) *dt * tmp%shell_vel ( :, shell_index )
tmp%core_vel (:,shell_index) = &
core_particle_set(shell_index)%v(:) * tmp%scale_v(:) * tmp%scale_v(:) +&
tmp%scale_v(:) * tmp%poly_v(:) *dmc * core_particle_set(shell_index)%f(:)
tmp%core_pos (:,shell_index) = &
core_particle_set(shell_index)%r(:) * tmp%scale_r(:) * tmp%scale_r(:) +&
tmp%scale_r(:) * tmp%poly_r(:) *dt * tmp%core_vel ( :, shell_index )
tmp%vel (:,iparticle) = fac_masss*tmp%shell_vel (:,shell_index)+&
fac_massc*tmp%core_vel (:,shell_index)
tmp%pos (:,iparticle) = fac_masss*tmp%shell_pos (:,shell_index)+&
fac_massc*tmp%core_pos (:,shell_index)
tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel(1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle)))
tmp%max_vel_sc = MAX(tmp%max_vel_sc,ABS(tmp%shell_vel(1,shell_index)-tmp%core_vel(1,shell_index)),&
ABS(tmp%shell_vel(2,shell_index)-tmp%core_vel(2,shell_index)),&
ABS(tmp%shell_vel(3,shell_index)-tmp%core_vel(3,shell_index)))
tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),&
ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),&
ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel,ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),&
ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),&
ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle)))
dsc(:) = tmp%shell_pos (:,shell_index)-tmp%core_pos (:,shell_index)-&
shell_particle_set(shell_index)%r(:)+core_particle_set(shell_index)%r(:)
tmp%max_dsc = MAX(tmp%max_dsc,ABS(dsc(1)),ABS(dsc(2)),ABS(dsc(3)))
dvsc(:) = tmp%shell_vel(:,shell_index)-tmp%core_vel(:,shell_index)- &
shell_particle_set(shell_index)%v(:)+core_particle_set(shell_index)%v(:)
tmp%max_dvel_sc = MAX(tmp%max_dvel_sc,ABS(dvsc(1)),ABS(dvsc(2)),ABS(dvsc(3)))
END DO ! iparticle_local
END IF
ELSE
nparticle_local = local_particles%n_el(iparticle_kind)
IF(PRESENT(u)) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
! Transform positions and velocities and forces
CALL transform_first(particle_set,tmp%pos,tmp%vel,&
iparticle,u,dm,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r)
tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel (1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle)))
tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),&
ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),&
ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel,ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),&
ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),&
ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle)))
END DO ! iparticle_local
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
tmp%vel (1,iparticle) = &
particle_set(iparticle)%v(1) * tmp%scale_v(1) * tmp%scale_v(1)+&
tmp%scale_v(1) * tmp%poly_v(1) * dm * particle_set(iparticle)%f(1)
tmp%vel (2,iparticle) = &
particle_set(iparticle)%v(2) * tmp%scale_v(2) * tmp%scale_v(2)+&
tmp%scale_v(2) * tmp%poly_v(2) * dm * particle_set(iparticle)%f(2)
tmp%vel (3,iparticle) = &
particle_set(iparticle)%v(3) * tmp%scale_v(3) * tmp%scale_v(3)+&
tmp%scale_v(3) * tmp%poly_v(3) * dm * particle_set(iparticle)%f(3)
tmp%pos (1,iparticle) = &
particle_set(iparticle)%r(1) * tmp%scale_r(1) * tmp%scale_r(1)+&
tmp%scale_r(1) * tmp%poly_r(1) * dt * tmp%vel ( 1, iparticle )
tmp%pos (2,iparticle) = &
particle_set(iparticle)%r(2) * tmp%scale_r(2) * tmp%scale_r(2)+&
tmp%scale_r(2) * tmp%poly_r(2) * dt * tmp%vel ( 2, iparticle )
tmp%pos (3,iparticle) = &
particle_set(iparticle)%r(3) * tmp%scale_r(3) * tmp%scale_r(3)+&
tmp%scale_r(3) * tmp%poly_r(3) * dt * tmp%vel ( 3, iparticle )
tmp%max_vel = MAX(tmp%max_vel,ABS(tmp%vel (1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle)))
tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),&
ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),&
ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel,ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),&
ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),&
ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle)))
END DO ! iparticle_local
nparticle_local = local_particles%n_el(iparticle_kind)
IF(PRESENT(u)) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
! Transform positions and velocities and forces
CALL transform_first(particle_set,tmp%pos,tmp%vel,&
iparticle,u,dm,dt,tmp%poly_v,tmp%poly_r,tmp%scale_v,tmp%scale_r)
tmp%max_vel = MAX(tmp%max_vel,&
ABS(tmp%vel (1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle)))
tmp%max_dr = MAX(tmp%max_dr,ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),&
ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),&
ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel,&
ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),&
ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),&
ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle)))
END DO ! iparticle_local
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
tmp%vel (1,iparticle) = &
particle_set(iparticle)%v(1) * tmp%scale_v(1) * tmp%scale_v(1)+&
tmp%scale_v(1) * tmp%poly_v(1) * dm * particle_set(iparticle)%f(1)
tmp%vel (2,iparticle) = &
particle_set(iparticle)%v(2) * tmp%scale_v(2) * tmp%scale_v(2)+&
tmp%scale_v(2) * tmp%poly_v(2) * dm * particle_set(iparticle)%f(2)
tmp%vel (3,iparticle) = &
particle_set(iparticle)%v(3) * tmp%scale_v(3) * tmp%scale_v(3)+&
tmp%scale_v(3) * tmp%poly_v(3) * dm * particle_set(iparticle)%f(3)
tmp%pos (1,iparticle) = &
particle_set(iparticle)%r(1) * tmp%scale_r(1) * tmp%scale_r(1)+&
tmp%scale_r(1) * tmp%poly_r(1) * dt * tmp%vel ( 1, iparticle )
tmp%pos (2,iparticle) = &
particle_set(iparticle)%r(2) * tmp%scale_r(2) * tmp%scale_r(2)+&
tmp%scale_r(2) * tmp%poly_r(2) * dt * tmp%vel ( 2, iparticle )
tmp%pos (3,iparticle) = &
particle_set(iparticle)%r(3) * tmp%scale_r(3) * tmp%scale_r(3)+&
tmp%scale_r(3) * tmp%poly_r(3) * dt * tmp%vel ( 3, iparticle )
tmp%max_vel = MAX(tmp%max_vel,&
ABS(tmp%vel (1,iparticle)),ABS(tmp%vel (2,iparticle)),ABS(tmp%vel (3,iparticle)))
tmp%max_dr = MAX(tmp%max_dr,&
ABS(particle_set(iparticle)%r(1)-tmp%pos (1,iparticle)),&
ABS(particle_set(iparticle)%r(2)-tmp%pos (2,iparticle)),&
ABS(particle_set(iparticle)%r(3)-tmp%pos (3,iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel,&
ABS(particle_set(iparticle)%v(1)-tmp%vel (1,iparticle)),&
ABS(particle_set(iparticle)%v(2)-tmp%vel (2,iparticle)),&
ABS(particle_set(iparticle)%v(3)-tmp%vel (3,iparticle)))
END DO
END IF
END IF
END IF
END DO ! iparticle_kind
END DO
END SUBROUTINE vv_first
! *****************************************************************************
@ -1191,90 +1208,91 @@ CONTAINS
atomic_kind => atomic_kind_set(iparticle_kind)
CALL get_atomic_kind(atomic_kind=atomic_kind,mass=mass,&
shell_active=is_shell)
dm = 0.5_dp*dt/mass
IF(is_shell .AND. shell_adiabatic) THEN
CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
dms = 0.5_dp*dt/shell%mass_shell
dmc = 0.5_dp*dt/shell%mass_core
fac_masss = shell%mass_shell/mass
fac_massc = shell%mass_core/mass
nparticle_local = local_particles%n_el(iparticle_kind)
IF(PRESENT(u)) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
! Transform velocities and forces shell
CALL transform_second(shell_particle_set,tmp%shell_vel,shell_index,&
u,dms,tmp%poly_v,tmp%scale_v)
IF (mass/=0.0_dp) THEN
dm = 0.5_dp*dt/mass
IF(is_shell .AND. shell_adiabatic) THEN
CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
dms = 0.5_dp*dt/shell%mass_shell
dmc = 0.5_dp*dt/shell%mass_core
fac_masss = shell%mass_shell/mass
fac_massc = shell%mass_core/mass
nparticle_local = local_particles%n_el(iparticle_kind)
IF(PRESENT(u)) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
! Transform velocities and forces shell
CALL transform_second(shell_particle_set,tmp%shell_vel,shell_index,&
u,dms,tmp%poly_v,tmp%scale_v)
! Transform velocities and forces core
CALL transform_second(core_particle_set,tmp%core_vel,shell_index,&
u,dmc,tmp%poly_v,tmp%scale_v)
! Transform velocities and forces core
CALL transform_second(core_particle_set,tmp%core_vel,shell_index,&
u,dmc,tmp%poly_v,tmp%scale_v)
! Derive velocties of the COM
tmp%vel(1,iparticle) = fac_masss*tmp%shell_vel(1,shell_index) + &
fac_massc*tmp%core_vel(1,shell_index)
tmp%vel(2,iparticle) = fac_masss*tmp%shell_vel(2,shell_index)+ &
fac_massc*tmp%core_vel(2,shell_index)
tmp%vel(3,iparticle) = fac_masss*tmp%shell_vel(3,shell_index) + &
fac_massc*tmp%core_vel(3,shell_index)
END DO ! iparticle_local
! Derive velocties of the COM
tmp%vel(1,iparticle) = fac_masss*tmp%shell_vel(1,shell_index) + &
fac_massc*tmp%core_vel(1,shell_index)
tmp%vel(2,iparticle) = fac_masss*tmp%shell_vel(2,shell_index)+ &
fac_massc*tmp%core_vel(2,shell_index)
tmp%vel(3,iparticle) = fac_masss*tmp%shell_vel(3,shell_index) + &
fac_massc*tmp%core_vel(3,shell_index)
END DO ! iparticle_local
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
tmp%shell_vel (1,shell_index) = &
tmp%shell_vel (1,shell_index) * tmp%scale_v(1) * tmp%scale_v(1) + &
tmp%scale_v(1) * tmp%poly_v(1) * dms*shell_particle_set(shell_index)%f(1)
tmp%shell_vel (2,shell_index) =&
tmp%shell_vel (2,shell_index) * tmp%scale_v(2) * tmp%scale_v(2) + &
tmp%scale_v(2) * tmp%poly_v(2) * dms*shell_particle_set(shell_index)%f(2)
tmp%shell_vel (3,shell_index) =&
tmp%shell_vel (3,shell_index) * tmp%scale_v(3) * tmp%scale_v(3) + &
tmp%scale_v(3) * tmp%poly_v(3) * dms*shell_particle_set(shell_index)%f(3)
tmp%core_vel (1,shell_index) =&
tmp%core_vel (1,shell_index) * tmp%scale_v(1) * tmp%scale_v(1) + &
tmp%scale_v(1) * tmp%poly_v(1) * dmc*core_particle_set(shell_index)%f(1)
tmp%core_vel (2,shell_index) =&
tmp%core_vel (2,shell_index) * tmp%scale_v(2) * tmp%scale_v(2)+ &
tmp%scale_v(2) * tmp%poly_v(2) * dmc*core_particle_set(shell_index)%f(2)
tmp%core_vel (3,shell_index) =&
tmp%core_vel (3,shell_index) * tmp%scale_v(3) * tmp%scale_v(3)+ &
tmp%scale_v(3) * tmp%poly_v(3) * dmc*core_particle_set(shell_index)%f(3)
tmp%vel ( 1, iparticle) = fac_masss*tmp%shell_vel (1,shell_index)+ &
fac_massc*tmp%core_vel (1,shell_index)
tmp%vel ( 2, iparticle) = fac_masss*tmp%shell_vel (2,shell_index)+ &
fac_massc*tmp%core_vel (2,shell_index)
tmp%vel ( 3, iparticle) = fac_masss*tmp%shell_vel (3,shell_index)+ &
fac_massc*tmp%core_vel (3,shell_index)
END DO ! iparticle_local
END IF
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
tmp%shell_vel (1,shell_index) = &
tmp%shell_vel (1,shell_index) * tmp%scale_v(1) * tmp%scale_v(1) + &
tmp%scale_v(1) * tmp%poly_v(1) * dms*shell_particle_set(shell_index)%f(1)
tmp%shell_vel (2,shell_index) =&
tmp%shell_vel (2,shell_index) * tmp%scale_v(2) * tmp%scale_v(2) + &
tmp%scale_v(2) * tmp%poly_v(2) * dms*shell_particle_set(shell_index)%f(2)
tmp%shell_vel (3,shell_index) =&
tmp%shell_vel (3,shell_index) * tmp%scale_v(3) * tmp%scale_v(3) + &
tmp%scale_v(3) * tmp%poly_v(3) * dms*shell_particle_set(shell_index)%f(3)
tmp%core_vel (1,shell_index) =&
tmp%core_vel (1,shell_index) * tmp%scale_v(1) * tmp%scale_v(1) + &
tmp%scale_v(1) * tmp%poly_v(1) * dmc*core_particle_set(shell_index)%f(1)
tmp%core_vel (2,shell_index) =&
tmp%core_vel (2,shell_index) * tmp%scale_v(2) * tmp%scale_v(2)+ &
tmp%scale_v(2) * tmp%poly_v(2) * dmc*core_particle_set(shell_index)%f(2)
tmp%core_vel (3,shell_index) =&
tmp%core_vel (3,shell_index) * tmp%scale_v(3) * tmp%scale_v(3)+ &
tmp%scale_v(3) * tmp%poly_v(3) * dmc*core_particle_set(shell_index)%f(3)
tmp%vel ( 1, iparticle) = fac_masss*tmp%shell_vel (1,shell_index)+ &
fac_massc*tmp%core_vel (1,shell_index)
tmp%vel ( 2, iparticle) = fac_masss*tmp%shell_vel (2,shell_index)+ &
fac_massc*tmp%core_vel (2,shell_index)
tmp%vel ( 3, iparticle) = fac_masss*tmp%shell_vel (3,shell_index)+ &
fac_massc*tmp%core_vel (3,shell_index)
END DO ! iparticle_local
END IF
ELSE
nparticle_local = local_particles%n_el(iparticle_kind)
IF(PRESENT(u)) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
CALL transform_second(particle_set,tmp%vel,iparticle,&
u,dm,tmp%poly_v,tmp%scale_v)
END DO ! iparticle_local
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
tmp%vel ( 1, iparticle) = &
tmp%vel ( 1, iparticle) * tmp%scale_v(1) * tmp%scale_v(1) + &
tmp%scale_v(1) * tmp%poly_v(1) * dm*particle_set(iparticle)%f(1)
tmp%vel ( 2, iparticle) = &
tmp%vel ( 2, iparticle) * tmp%scale_v(2) * tmp%scale_v(2) + &
tmp%scale_v(2) * tmp%poly_v(2) * dm*particle_set(iparticle)%f(2)
tmp%vel ( 3, iparticle) = &
tmp%vel ( 3, iparticle) * tmp%scale_v(3) * tmp%scale_v(3) + &
tmp%scale_v(3) * tmp%poly_v(3) * dm*particle_set(iparticle)%f(3)
END DO ! iparticle_local
nparticle_local = local_particles%n_el(iparticle_kind)
IF(PRESENT(u)) THEN
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
CALL transform_second(particle_set,tmp%vel,iparticle,&
u,dm,tmp%poly_v,tmp%scale_v)
END DO ! iparticle_local
ELSE
DO iparticle_local=1,nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
tmp%vel ( 1, iparticle) = &
tmp%vel ( 1, iparticle) * tmp%scale_v(1) * tmp%scale_v(1) + &
tmp%scale_v(1) * tmp%poly_v(1) * dm*particle_set(iparticle)%f(1)
tmp%vel ( 2, iparticle) = &
tmp%vel ( 2, iparticle) * tmp%scale_v(2) * tmp%scale_v(2) + &
tmp%scale_v(2) * tmp%poly_v(2) * dm*particle_set(iparticle)%f(2)
tmp%vel ( 3, iparticle) = &
tmp%vel ( 3, iparticle) * tmp%scale_v(3) * tmp%scale_v(3) + &
tmp%scale_v(3) * tmp%poly_v(3) * dm*particle_set(iparticle)%f(3)
END DO
END IF
END IF
END IF
END DO ! iparticle_kind
END DO
END SUBROUTINE vv_second

62
src/md_vel_utils.F
View file

@ -9,6 +9,8 @@
!> \author CJM
!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008
!> reorganization of the original routines/modules
!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
!> (patch by Marcel Baer)
! *****************************************************************************
MODULE md_vel_utils
USE atomic_kind_list_types, ONLY: atomic_kind_list_type
@ -142,13 +144,15 @@ CONTAINS
DO i = 1, SIZE(part)
atomic_kind => part(i)%atomic_kind
CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
SELECT CASE(is_fixed(i))
CASE(use_perd_x,use_perd_y,use_perd_z,use_perd_xy,use_perd_xz,use_perd_yz,use_perd_none)
vcom(1) = vcom(1) + part(i)%v(1) * mass
vcom(2) = vcom(2) + part(i)%v(2) * mass
vcom(3) = vcom(3) + part(i)%v(3) * mass
denom = denom + mass
END SELECT
IF (mass.NE.0.0) THEN
SELECT CASE(is_fixed(i))
CASE(use_perd_x,use_perd_y,use_perd_z,use_perd_xy,use_perd_xz,use_perd_yz,use_perd_none)
vcom(1) = vcom(1) + part(i)%v(1) * mass
vcom(2) = vcom(2) + part(i)%v(2) * mass
vcom(3) = vcom(3) + part(i)%v(3) * mass
denom = denom + mass
END SELECT
END IF
END DO
vcom = vcom/denom
IF (PRESENT(ecom)) THEN
@ -707,27 +711,29 @@ CONTAINS
part(i)%v(1) = 0.0_dp
part(i)%v(2) = 0.0_dp
part(i)%v(3) = 0.0_dp
SELECT CASE(is_fixed(i))
CASE (use_perd_x)
part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_y)
part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_z)
part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_xy)
part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_xz)
part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_yz)
part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_none)
part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
END SELECT
IF (mass.NE.0.0) THEN
SELECT CASE(is_fixed(i))
CASE (use_perd_x)
part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_y)
part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_z)
part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_xy)
part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_xz)
part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_yz)
part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
CASE (use_perd_none)
part(i)%v(1) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(2) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
part(i)%v(3) = next_random_number(globenv%gaussian_rng_stream,error=error) / SQRT(mass)
END SELECT
END IF
END DO
! Subtract the vcom

60
src/molecule_kind_types.F
View file

@ -7,7 +7,8 @@
!> \brief Define the molecule kind structure types and the corresponding
!> functionality
!> \par History
!> none
!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
!> (patch by Marcel Baer)
!> \author MK (22.08.2003)
! *****************************************************************************
MODULE molecule_kind_types
@ -138,6 +139,13 @@ MODULE molecule_kind_types
TYPE(restraint_type) :: restraint
END TYPE g4x6_constraint_type
! *****************************************************************************
TYPE vsite_constraint_type
INTEGER :: a,b,c,d
REAL(KIND = dp) :: wbc, wdc
TYPE(restraint_type) :: restraint
END TYPE vsite_constraint_type
! *****************************************************************************
TYPE fixd_constraint_type
TYPE(restraint_type) :: restraint
@ -166,9 +174,10 @@ MODULE molecule_kind_types
TYPE(onfo_kind_type), DIMENSION(:), POINTER :: onfo_kind_set
TYPE(onfo_type), DIMENSION(:), POINTER :: onfo_list
TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list
TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list
TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
TYPE(shell_type), DIMENSION(:), POINTER :: shell_list
CHARACTER(LEN=default_string_length) :: name
REAL(KIND = dp) :: charge,&
@ -183,6 +192,7 @@ MODULE molecule_kind_types
nub,&
ng3x3, ng3x3_restraint,&
ng4x6, ng4x6_restraint,&
nvsite, nvsite_restraint,&
nfixd, nfixd_restraint,&
nmolecule, nshell
TYPE(colvar_counters) :: ncolv
@ -214,6 +224,7 @@ MODULE molecule_kind_types
colvar_constraint_type,&
g3x3_constraint_type,&
g4x6_constraint_type,&
vsite_constraint_type,&
fixd_constraint_type,&
local_fixd_constraint_type,&
molecule_kind_type,&
@ -327,6 +338,7 @@ CONTAINS
NULLIFY (molecule_kind_set(imolecule_kind)%onfo_kind_set)
NULLIFY (molecule_kind_set(imolecule_kind)%g3x3_list)
NULLIFY (molecule_kind_set(imolecule_kind)%g4x6_list)
NULLIFY (molecule_kind_set(imolecule_kind)%vsite_list)
NULLIFY (molecule_kind_set(imolecule_kind)%fixd_list)
NULLIFY (molecule_kind_set(imolecule_kind)%shell_list)
NULLIFY (molecule_kind_set(imolecule_kind)%torsion_list)
@ -348,9 +360,11 @@ CONTAINS
molecule_kind_set(imolecule_kind)%ncolv,error)
molecule_kind_set(imolecule_kind)%ng3x3 = 0
molecule_kind_set(imolecule_kind)%ng4x6 = 0
molecule_kind_set(imolecule_kind)%nvsite = 0
molecule_kind_set(imolecule_kind)%nfixd = 0
molecule_kind_set(imolecule_kind)%ng3x3_restraint = 0
molecule_kind_set(imolecule_kind)%ng4x6_restraint = 0
molecule_kind_set(imolecule_kind)%nvsite_restraint = 0
molecule_kind_set(imolecule_kind)%nfixd_restraint = 0
molecule_kind_set(imolecule_kind)%nmolecule = 0
molecule_kind_set(imolecule_kind)%ntorsion = 0
@ -447,6 +461,10 @@ CONTAINS
DEALLOCATE (molecule_kind_set(imolecule_kind)%g4x6_list,STAT=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%vsite_list)) THEN
DEALLOCATE (molecule_kind_set(imolecule_kind)%vsite_list,STAT=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%fixd_list)) THEN
DEALLOCATE (molecule_kind_set(imolecule_kind)%fixd_list,STAT=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
@ -489,11 +507,11 @@ CONTAINS
! *****************************************************************************
SUBROUTINE get_molecule_kind(molecule_kind,atom_list,bond_list,bend_list,&
ub_list,impr_list,onfo_list,colv_list,fixd_list,&
g3x3_list,g4x6_list,torsion_list,shell_list,&
g3x3_list,g4x6_list,vsite_list,torsion_list,shell_list,&
name,mass,charge,kind_number,natom,nbend,nbond,nub,&
nimpr,nonfo,nconstraint,nconstraint_fixd,nfixd,ncolv,ng3x3,ng4x6,&
nfixd_restraint,ng3x3_restraint, ng4x6_restraint,&
nrestraints,nmolecule,nsgf,nshell,ntorsion,&
nvsite,nfixd_restraint,ng3x3_restraint, ng4x6_restraint,&
nvsite_restraint,nrestraints,nmolecule,nsgf,nshell,ntorsion,&
molecule_list,nelectron, bond_kind_set,bend_kind_set,&
ub_kind_set,impr_kind_set,onfo_kind_set,torsion_kind_set,&
molname_generated)
@ -519,6 +537,8 @@ CONTAINS
DIMENSION(:), OPTIONAL, POINTER :: g3x3_list
TYPE(g4x6_constraint_type), &
DIMENSION(:), OPTIONAL, POINTER :: g4x6_list
TYPE(vsite_constraint_type), &
DIMENSION(:), OPTIONAL, POINTER :: vsite_list
TYPE(torsion_type), DIMENSION(:), &
OPTIONAL, POINTER :: torsion_list
TYPE(shell_type), DIMENSION(:), &
@ -532,9 +552,9 @@ CONTAINS
nconstraint_fixd, nfixd
TYPE(colvar_counters), INTENT(out), &
OPTIONAL :: ncolv
INTEGER, INTENT(OUT), OPTIONAL :: ng3x3, ng4x6, nfixd_restraint, &
ng3x3_restraint, ng4x6_restraint, nrestraints, nmolecule, nsgf, nshell, &
ntorsion
INTEGER, INTENT(OUT), OPTIONAL :: ng3x3, ng4x6, nvsite, nfixd_restraint, &
ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, &
nmolecule, nsgf, nshell, ntorsion
INTEGER, DIMENSION(:), OPTIONAL, POINTER :: molecule_list
INTEGER, INTENT(OUT), OPTIONAL :: nelectron
TYPE(bond_kind_type), DIMENSION(:), &
@ -573,6 +593,7 @@ CONTAINS
IF (PRESENT(colv_list)) colv_list => molecule_kind%colv_list
IF (PRESENT(g3x3_list)) g3x3_list => molecule_kind%g3x3_list
IF (PRESENT(g4x6_list)) g4x6_list => molecule_kind%g4x6_list
IF (PRESENT(vsite_list)) vsite_list => molecule_kind%vsite_list
IF (PRESENT(fixd_list)) fixd_list => molecule_kind%fixd_list
IF (PRESENT(torsion_list)) torsion_list => molecule_kind%torsion_list
IF (PRESENT(shell_list)) shell_list => molecule_kind%shell_list
@ -589,10 +610,12 @@ CONTAINS
IF (PRESENT(nonfo)) nonfo = molecule_kind%nonfo
IF (PRESENT(nconstraint)) nconstraint = (molecule_kind%ncolv%ntot - molecule_kind%ncolv%nrestraint) +&
3*(molecule_kind%ng3x3-molecule_kind%ng3x3_restraint) +&
6*(molecule_kind%ng4x6-molecule_kind%ng4x6_restraint)
6*(molecule_kind%ng4x6-molecule_kind%ng4x6_restraint)+&
3*(molecule_kind%nvsite-molecule_kind%nvsite_restraint)
IF (PRESENT(ncolv)) ncolv = molecule_kind%ncolv
IF (PRESENT(ng3x3)) ng3x3= molecule_kind%ng3x3
IF (PRESENT(ng4x6)) ng4x6= molecule_kind%ng4x6
IF (PRESENT(nvsite)) nvsite= molecule_kind%nvsite
! Number of atoms that have one or more components fixed
IF (PRESENT(nfixd)) nfixd= molecule_kind%nfixd
! Number of degrees of freedom fixed
@ -614,10 +637,12 @@ CONTAINS
END IF
IF (PRESENT(ng3x3_restraint)) ng3x3_restraint= molecule_kind%ng3x3_restraint
IF (PRESENT(ng4x6_restraint)) ng4x6_restraint= molecule_kind%ng4x6_restraint
IF (PRESENT(nvsite_restraint)) nvsite_restraint= molecule_kind%nvsite_restraint
IF (PRESENT(nfixd_restraint)) nfixd_restraint= molecule_kind%nfixd_restraint
IF (PRESENT(nrestraints)) nrestraints = molecule_kind%ncolv%nrestraint +&
molecule_kind%ng3x3_restraint +&
molecule_kind%ng4x6_restraint
molecule_kind%ng4x6_restraint+&
molecule_kind%nvsite_restraint
IF (PRESENT(nmolecule)) nmolecule = molecule_kind%nmolecule
IF (PRESENT(nshell)) nshell= molecule_kind%nshell
IF (PRESENT(ntorsion)) ntorsion= molecule_kind%ntorsion
@ -802,8 +827,8 @@ CONTAINS
nbend,bend_list,nub,ub_list,nimpr,impr_list,&
nonfo,onfo_list,ntorsion,torsion_list,fixd_list,&
ncolv,colv_list,ng3x3,g3x3_list,ng4x6,nfixd,g4x6_list,&
ng3x3_restraint,ng4x6_restraint,nfixd_restraint,&
nshell, shell_list,&
nvsite,vsite_list,ng3x3_restraint,ng4x6_restraint,nfixd_restraint,&
nshell, shell_list,nvsite_restraint,&
bond_kind_set,bend_kind_set,ub_kind_set,&
torsion_kind_set,impr_kind_set,onfo_kind_set,&
nelectron,nsgf,molname_generated)
@ -845,11 +870,15 @@ CONTAINS
INTEGER, INTENT(IN), OPTIONAL :: ng4x6, nfixd
TYPE(g4x6_constraint_type), &
DIMENSION(:), OPTIONAL, POINTER :: g4x6_list
INTEGER, INTENT(IN), OPTIONAL :: nvsite
TYPE(vsite_constraint_type), &
DIMENSION(:), OPTIONAL, POINTER :: vsite_list
INTEGER, INTENT(IN), OPTIONAL :: ng3x3_restraint, &
ng4x6_restraint, &
nfixd_restraint, nshell
TYPE(shell_type), DIMENSION(:), &
OPTIONAL, POINTER :: shell_list
INTEGER, INTENT(IN), OPTIONAL :: nvsite_restraint
TYPE(bond_kind_type), DIMENSION(:), &
OPTIONAL, POINTER :: bond_kind_set
TYPE(bend_kind_type), DIMENSION(:), &
@ -901,11 +930,14 @@ CONTAINS
IF (PRESENT(ng3x3)) molecule_kind%ng3x3 = ng3x3
IF (PRESENT(g3x3_list)) molecule_kind%g3x3_list => g3x3_list
IF (PRESENT(ng4x6)) molecule_kind%ng4x6 = ng4x6
IF (PRESENT(nvsite)) molecule_kind%nvsite = nvsite
IF (PRESENT(nfixd)) molecule_kind%nfixd = nfixd
IF (PRESENT(nfixd_restraint)) molecule_kind%nfixd_restraint = nfixd_restraint
IF (PRESENT(ng3x3_restraint)) molecule_kind%ng3x3_restraint = ng3x3_restraint
IF (PRESENT(ng4x6_restraint)) molecule_kind%ng4x6_restraint = ng4x6_restraint
IF (PRESENT(nvsite_restraint)) molecule_kind%nvsite_restraint = nvsite_restraint
IF (PRESENT(g4x6_list)) molecule_kind%g4x6_list => g4x6_list
IF (PRESENT(vsite_list)) molecule_kind%vsite_list => vsite_list
IF (PRESENT(fixd_list)) molecule_kind%fixd_list => fixd_list
IF (PRESENT(bond_kind_set)) molecule_kind%bond_kind_set => bond_kind_set
IF (PRESENT(bend_kind_set)) molecule_kind%bend_kind_set => bend_kind_set

7
src/molecule_types_new.F
View file

@ -7,6 +7,8 @@
!> \brief Define the data structure for the molecule information.
!> \par History
!> JGH (22.05.2004) add last_atom information
!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
!> (patch by Marcel Baer)
!> \author MK (29.08.2003)
! *****************************************************************************
MODULE molecule_types_new
@ -21,7 +23,8 @@ MODULE molecule_types_new
fixd_constraint_type,&
g3x3_constraint_type,&
g4x6_constraint_type,&
molecule_kind_type
molecule_kind_type,&
vsite_constraint_type
USE termination, ONLY: stop_memory,&
stop_program
#include "cp_common_uses.h"
@ -86,10 +89,12 @@ MODULE molecule_types_new
INTEGER :: ntot, nrestraint
INTEGER :: ng3x3, ng3x3_restraint
INTEGER :: ng4x6, ng4x6_restraint
INTEGER :: nvsite, nvsite_restraint
TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list
TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list
TYPE(local_colvar_constraint_type), DIMENSION ( : ) , POINTER :: lcolv
TYPE(local_g3x3_constraint_type), DIMENSION ( : ) , POINTER :: lg3x3
TYPE(local_g4x6_constraint_type), DIMENSION ( : ) , POINTER :: lg4x6

100
src/topology_constraint_util.F
View file

@ -35,7 +35,8 @@ MODULE topology_constraint_util
USE molecule_kind_types, ONLY: &
atom_type, bond_type, colvar_constraint_type, fixd_constraint_type, &
g3x3_constraint_type, g4x6_constraint_type, get_molecule_kind, &
molecule_kind_type, set_molecule_kind, setup_colvar_counters
molecule_kind_type, set_molecule_kind, setup_colvar_counters, &
vsite_constraint_type
USE molecule_types_new, ONLY: get_molecule,&
global_constraint_type,&
local_colvar_constraint_type,&
@ -97,7 +98,7 @@ CONTAINS
k1loc, k2loc, kk, last, last_atom, m, n_start_colv, natom, nbond, &
ncolv_glob, ncolv_mol, nfixd_list_gci, nfixd_restart, nfixd_restraint, &
nfixed_atoms, ng3x3, ng3x3_restraint, ng4x6, ng4x6_restraint, nhdist, &
nmolecule, nrep, offset, output_unit, stat
nmolecule, nrep, nvsite, nvsite_restraint, offset, output_unit, stat
INTEGER, DIMENSION(:), POINTER :: constr_x_glob, inds, &
molecule_list
LOGICAL :: exclude_mm, exclude_qm, failure, fix_atom_mm, &
@ -134,6 +135,8 @@ CONTAINS
TYPE(molecule_type), POINTER :: molecule
TYPE(section_vals_type), POINTER :: colvar_rest, fixd_restr_rest, &
hbonds_section
TYPE(vsite_constraint_type), &
DIMENSION(:), POINTER :: vsite_list
failure = .FALSE.
NULLIFY(logger, constr_x_mol, constr_x_glob)
@ -179,15 +182,19 @@ CONTAINS
gci%fixd_list,&
gci%colv_list,&
gci%g3x3_list,&
gci%g4x6_list)
gci%g4x6_list,&
gci%vsite_list)
gci%ntot = 0
gci%ng3x3 = 0
gci%ng4x6 = 0
gci%nvsite = 0
gci%ng3x3_restraint = 0
gci%ng4x6_restraint = 0
gci%nvsite_restraint = 0
CALL setup_colvar_counters(gci%colv_list,gci%ncolv,error)
gci%nrestraint = gci%ng3x3_restraint +&
gci%ng4x6_restraint +&
gci%nvsite_restraint +&
gci%ncolv%nrestraint
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
@ -574,7 +581,54 @@ CONTAINS
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! 6. Set the group fixed_atom constraint fixd_list
! 6. Set the group vsite constraint vsite_list
!-----------------------------------------------------------------------------
IF(topology%const_vsite) THEN
CALL give_constraint_array(cons_info%const_vsite_mol,&
cons_info%const_vsite_molname,&
cons_info%vsite_intermolecular,&
constr_x_mol,&
constr_x_glob,&
molecule_kind_set,&
cons_info%vsite_exclude_qm,&
cons_info%vsite_exclude_mm,&
error)
! Intramolecular constraints
DO ii = 1, SIZE(molecule_kind_set)
molecule_kind => molecule_kind_set(ii)
CALL get_molecule_kind(molecule_kind=molecule_kind,&
nmolecule=nmolecule,molecule_list=molecule_list)
nvsite = SIZE(constr_x_mol(ii)%constr)
ALLOCATE(vsite_list(nvsite),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
CALL setup_vsite_list(vsite_list, constr_x_mol(ii)%constr, cons_info, nvsite_restraint, error)
CALL set_molecule_kind(molecule_kind,nvsite=nvsite,nvsite_restraint=nvsite_restraint,&
vsite_list=vsite_list)
END DO
DO j = 1, SIZE(constr_x_mol)
DEALLOCATE(constr_x_mol(j)%constr,stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END DO
DEALLOCATE(constr_x_mol,stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
! Intermolecular constraints
IF (ASSOCIATED(constr_x_glob)) THEN
nvsite = SIZE(constr_x_glob)
ALLOCATE(vsite_list(nvsite),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
CALL setup_vsite_list(vsite_list, constr_x_glob, cons_info, nvsite_restraint, error)
gci%vsite_list => vsite_list
gci%nvsite = nvsite
gci%nvsite_restraint = nvsite_restraint
! Total number of Intermolecular constraints
gci%ntot = gci%nvsite + gci%ntot
DEALLOCATE(constr_x_glob,stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
END IF
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! 7. Set the group fixed_atom constraint fixd_list
!-----------------------------------------------------------------------------
IF(topology%const_atom) THEN
ALLOCATE(fixd_list_gci(SIZE(particle_set)),stat=stat)
@ -846,6 +900,7 @@ CONTAINS
! Final setup of the number of possible restraints
gci%nrestraint = gci%ng3x3_restraint +&
gci%ng4x6_restraint +&
gci%nvsite_restraint +&
gci%ncolv%nrestraint
CALL cp_print_key_finished_output(iw,logger,subsys_section,&
"PRINT%TOPOLOGY_INFO/UTIL_INFO",error=error)
@ -1013,6 +1068,43 @@ CONTAINS
END SUBROUTINE setup_g4x6_list
! *****************************************************************************
!> \brief Setup the vsite_list for the packing of constraints
!> \par History
!> \author Marcel Baer [2008]
! *****************************************************************************
SUBROUTINE setup_vsite_list(vsite_list, ilist, cons_info, nvsite_restraint, error)
TYPE(vsite_constraint_type), &
DIMENSION(:), POINTER :: vsite_list
INTEGER, DIMENSION(:), POINTER :: ilist
TYPE(constraint_info_type), POINTER :: cons_info
INTEGER, INTENT(OUT) :: nvsite_restraint
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'setup_vsite_list', &
routineP = moduleN//':'//routineN
INTEGER :: j, nvsite
LOGICAL :: failure
failure = .FALSE.
nvsite = 0
nvsite_restraint = 0
DO nvsite = 1, SIZE(ilist)
j = ilist(nvsite)
vsite_list(nvsite)%a = cons_info%const_vsite_a(j)
vsite_list(nvsite)%b = cons_info%const_vsite_b(j)
vsite_list(nvsite)%c = cons_info%const_vsite_c(j)
vsite_list(nvsite)%d = cons_info%const_vsite_d(j)
vsite_list(nvsite)%wbc = cons_info%const_vsite_wbc(j)
vsite_list(nvsite)%wdc = cons_info%const_vsite_wdc(j)
! Restraint
vsite_list(nvsite)%restraint%active = cons_info%vsite_restraint(j)
vsite_list(nvsite)%restraint%k0 = cons_info%vsite_k0(j)
IF (vsite_list(nvsite)%restraint%active) nvsite_restraint = nvsite_restraint+ 1
END DO
END SUBROUTINE setup_vsite_list
! *****************************************************************************
!> \brief Setup the lcolv for the packing of constraints
!> \par History
!> Updated 2007 for intermolecular constraints

94
src/topology_input.F
View file

@ -126,10 +126,8 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro
LOGICAL :: explicit, failure
REAL(KIND=dp), DIMENSION(:), POINTER :: rlist
TYPE(constraint_info_type), POINTER :: cons_info
TYPE(section_vals_type), POINTER :: collective_section, &
fix_atom_section, &
g3x3_section, g4x6_section, &
hbonds_section
TYPE(section_vals_type), POINTER :: collective_section, fix_atom_section, &
g3x3_section, g4x6_section, hbonds_section, vsite_section
cons_info => topology%cons_info
failure = .FALSE.
@ -137,6 +135,7 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro
hbonds_section => section_vals_get_subs_vals(constraint_section,"HBONDS",error=error)
g3x3_section => section_vals_get_subs_vals(constraint_section,"G3X3",error=error)
g4x6_section => section_vals_get_subs_vals(constraint_section,"G4X6",error=error)
vsite_section => section_vals_get_subs_vals(constraint_section,"VIRTUAL_SITE",error=error)
fix_atom_section => section_vals_get_subs_vals(constraint_section,"FIXED_ATOMS",error=error)
collective_section => section_vals_get_subs_vals(constraint_section,"COLLECTIVE",error=error)
! HBONDS
@ -333,6 +332,93 @@ SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section, erro
cons_info%const_g46_dcd(ig) = rlist(6)
END DO
END IF
! virtual
CALL section_vals_get(vsite_section,explicit=explicit, n_repetition=ncons, error=error)
IF (explicit) THEN
topology%const_vsite =.TRUE.
cons_info%nconst_vsite = ncons
!
ALLOCATE(cons_info%const_vsite_mol(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%const_vsite_molname(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%const_vsite_a(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%const_vsite_b(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%const_vsite_c(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%const_vsite_d(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%const_vsite_wbc(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%const_vsite_wdc(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%vsite_intermolecular(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%vsite_restraint(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%vsite_k0(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%vsite_exclude_qm(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(cons_info%vsite_exclude_mm(ncons),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
DO ig = 1, ncons
CALL check_restraint(vsite_section,&
is_restraint=cons_info%vsite_restraint(ig), &
k0=cons_info%vsite_k0(ig),&
i_rep_section=ig,&
label="Virtual_SITE",&
error=error)
cons_info%const_vsite_mol(ig) = 0
cons_info%const_vsite_molname(ig) = "UNDEF"
! Exclude QM or MM
CALL section_vals_val_get(vsite_section,"EXCLUDE_QM",i_rep_section=ig,&
l_val=cons_info%vsite_exclude_qm(ig),error=error)
CALL section_vals_val_get(vsite_section,"EXCLUDE_MM",i_rep_section=ig,&
l_val=cons_info%vsite_exclude_mm(ig),error=error)
! Intramolecular restraint
CALL section_vals_val_get(vsite_section,"INTERMOLECULAR",i_rep_section=ig,&
l_val=cons_info%vsite_intermolecular(ig),error=error)
! If it is intramolecular let's unset (in case user did it)
! the molecule and molname field
IF (cons_info%vsite_intermolecular(ig)) THEN
CALL section_vals_val_unset(vsite_section,"MOLECULE",i_rep_section=ig,error=error)
CALL section_vals_val_unset(vsite_section,"MOLNAME",i_rep_section=ig,error=error)
END IF
! Let's tag to which molecule we want to apply constraints
CALL section_vals_val_get(vsite_section,"MOLECULE",i_rep_section=ig,&
n_rep_val=nrep,error=error)
IF (nrep/=0) THEN
CALL section_vals_val_get(vsite_section,"MOLECULE",i_rep_section=ig,&
i_val=cons_info%const_vsite_mol(ig),error=error)
END IF
CALL section_vals_val_get(vsite_section,"MOLNAME",i_rep_section=ig,&
n_rep_val=nrep,error=error)
IF (nrep/=0) THEN
CALL section_vals_val_get(vsite_section,"MOLNAME",i_rep_section=ig,&
c_val=cons_info%const_vsite_molname(ig),error=error)
END IF
IF ((cons_info%const_vsite_mol(ig)/=0).AND.(cons_info%const_vsite_molname(ig)/="UNDEF")) THEN
CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure)
END IF
IF ((cons_info%const_vsite_mol(ig)==0).AND.(cons_info%const_vsite_molname(ig)=="UNDEF").AND.&
(.NOT.cons_info%vsite_intermolecular(ig))) THEN
CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure)
END IF
CALL section_vals_val_get(vsite_section,"ATOMS",i_rep_section=ig,&
i_vals=ilist,error=error)
CALL section_vals_val_get(vsite_section,"PARAMETERS",i_rep_section=ig,&
r_vals=rlist,error=error)
cons_info%const_vsite_a(ig) = ilist(1)
cons_info%const_vsite_b(ig) = ilist(2)
cons_info%const_vsite_c(ig) = ilist(3)
cons_info%const_vsite_d(ig) = ilist(4)
cons_info%const_vsite_wbc(ig) = rlist(1)
cons_info%const_vsite_wdc(ig) = rlist(2)
END DO
END IF
! FIXED ATOMS
CALL section_vals_get(fix_atom_section,explicit=explicit, n_repetition=ncons, error=error)
IF (explicit) THEN

85
src/topology_types.F
View file

@ -9,6 +9,8 @@
!> change name from input_file_name to coord_...
!> added topology file
!> added atom_names
!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
!> (patch by Marcel Baer)
!> \author CJM & JGH
! *****************************************************************************
MODULE topology_types
@ -119,6 +121,20 @@ MODULE topology_types
LOGICAL, POINTER :: g46_restraint(:) ! Restraints control
REAL(KIND=dp), POINTER :: g46_k0(:) ! Restraints control
LOGICAL, POINTER, DIMENSION(:) :: g46_exclude_qm, g46_exclude_mm
! virtual_site
INTEGER :: nconst_vsite
INTEGER, POINTER :: const_vsite_mol(:)
CHARACTER(LEN=default_string_length), POINTER :: const_vsite_molname(:)
INTEGER, POINTER :: const_vsite_a(:)
INTEGER, POINTER :: const_vsite_b(:)
INTEGER, POINTER :: const_vsite_c(:)
INTEGER, POINTER :: const_vsite_d(:)
REAL(KIND=dp), POINTER :: const_vsite_wbc(:)
REAL(KIND=dp), POINTER :: const_vsite_wdc(:)
LOGICAL, POINTER :: vsite_intermolecular(:)
LOGICAL, POINTER :: vsite_restraint(:) ! Restraints control
REAL(KIND=dp), POINTER :: vsite_k0(:) ! Restraints control
LOGICAL, POINTER, DIMENSION(:) :: vsite_exclude_qm, vsite_exclude_mm
END TYPE constraint_info_type
! *****************************************************************************
@ -151,6 +167,7 @@ MODULE topology_types
LOGICAL :: const_colv
LOGICAL :: const_33
LOGICAL :: const_46
LOGICAL :: const_vsite
LOGICAL :: charge_occup
LOGICAL :: charge_beta
LOGICAL :: charge_extended
@ -226,6 +243,7 @@ CONTAINS
topology%const_colv=.FALSE.
topology%const_33=.FALSE.
topology%const_46=.FALSE.
topology%const_vsite=.FALSE.
topology%charge_occup=.FALSE.
topology%charge_beta=.FALSE.
topology%charge_extended=.FALSE.
@ -360,6 +378,20 @@ CONTAINS
NULLIFY(constraint_info%g46_restraint)
NULLIFY(constraint_info%g46_k0)
NULLIFY(constraint_info%g46_exclude_qm, constraint_info%g46_exclude_mm)
! virtual_site
constraint_info%nconst_vsite = 0
NULLIFY(constraint_info%const_vsite_mol)
NULLIFY(constraint_info%const_vsite_molname)
NULLIFY(constraint_info%const_vsite_a)
NULLIFY(constraint_info%const_vsite_b)
NULLIFY(constraint_info%const_vsite_c)
NULLIFY(constraint_info%const_vsite_d)
NULLIFY(constraint_info%const_vsite_wbc)
NULLIFY(constraint_info%const_vsite_wdc)
NULLIFY(constraint_info%vsite_intermolecular)
NULLIFY(constraint_info%vsite_restraint)
NULLIFY(constraint_info%vsite_k0)
NULLIFY(constraint_info%vsite_exclude_qm, constraint_info%vsite_exclude_mm)
END SUBROUTINE init_constraint
@ -798,6 +830,59 @@ CONTAINS
DEALLOCATE(constraint_info%g46_exclude_mm, stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
! virtual_site
IF(ASSOCIATED(constraint_info%const_vsite_mol)) THEN
DEALLOCATE(constraint_info%const_vsite_mol,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%const_vsite_molname)) THEN
DEALLOCATE(constraint_info%const_vsite_molname,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%const_vsite_a)) THEN
DEALLOCATE(constraint_info%const_vsite_a,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%const_vsite_b)) THEN
DEALLOCATE(constraint_info%const_vsite_b,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%const_vsite_c)) THEN
DEALLOCATE(constraint_info%const_vsite_c,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%const_vsite_d)) THEN
DEALLOCATE(constraint_info%const_vsite_d,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%const_vsite_wbc)) THEN
DEALLOCATE(constraint_info%const_vsite_wbc,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%const_vsite_wdc)) THEN
DEALLOCATE(constraint_info%const_vsite_wdc,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%vsite_intermolecular)) THEN
DEALLOCATE(constraint_info%vsite_intermolecular,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%vsite_restraint)) THEN
DEALLOCATE(constraint_info%vsite_restraint,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%vsite_k0)) THEN
DEALLOCATE(constraint_info%vsite_k0,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%vsite_exclude_qm)) THEN
DEALLOCATE(constraint_info%vsite_exclude_qm, stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(constraint_info%vsite_exclude_mm)) THEN
DEALLOCATE(constraint_info%vsite_exclude_mm, stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
END IF
END SUBROUTINE deallocate_constraint

4
tests/Fist/regtest-10/TEST_FILES_RESET
View file

@ -31,3 +31,7 @@ nh3-meta-2.inp
# numerical instabilities
nh3-meta-3.inp
nh3-meta-4.inp
# numerics
nh3-meta-2.inp
nh3-meta-3.inp
nh3-meta-4.inp

5
tests/Fist/regtest-4/TEST_FILES_RESET
View file

@ -177,3 +177,8 @@ nh3-meta-walks_2.inp
nh3-meta-walks_1r.inp
# unified cube_center
nh3-meta-walks_2r.inp
# numerics
nh3-meta-walks_1.inp
nh3-meta-walks_1r.inp
nh3-meta-walks_2.inp
nh3-meta-walks_2r.inp

2
tests/Fist/regtest-5/TEST_FILES
View file

@ -40,3 +40,5 @@ H2OFe-5_75.inp 2
acn-2.inp 0
H2OFe-5_75-2.inp 0
I3.inp 2
#TIP4P
water_tip4p.inp 2

224
tests/Fist/regtest-5/water_tip4p.inp Normal file
View file

@ -0,0 +1,224 @@
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
&SPLINE
EMAX_SPLINE 1000000.0
RCUT_NB 12.0
&END SPLINE
!GROMACS parms
![ atomtypes ]
!; full atom descriptions are available in ffoplsaa.atp
!; name bond_type mass charge ptype sigma epsilon
! opls_113 OW 15.99940 0.000 A 3.15365e-01 6.48520e-01
! opls_114 HW 1.00800 0.520 A 0.00000e+00 0.00000e+00
! opls_115 MW 0.00000 -1.040 D 0.00000e+00 0.00000e+00
!
![ bondtypes ]
!; i j func b0 kb
! OW HW 1 0.09572 502080.0 ; For TIP4F Water - wlj 1/98
! OW LP 1 0.01750 753120.0 ; -idem-
!
![ angletypes ]
!; i j k func th0 cth
! HW OW HW 1 109.500 627.600 ; For TIP4F Water - wj 1/98
! HW OW LP 1 54.750 418.400 ; For TIP4F Water - wj 1/98
&BEND
ATOMS HW OW HW
THETA0 [deg] 109.500
!K [rad^2kjmol] 627.600
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS HW OW MW
THETA0 [deg] 54.750
!K [rad^2kjmol] 418.400
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS MW HW OW
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS MW HW HW
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS OW HW HW
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS OW MW HW
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BEND
ATOMS HW MW HW
THETA0 [deg] 109.500
K [rad^2kjmol] 0.000
KIND G87
&END BEND
&BOND
ATOMS OW HW
K [nm^-2kjmol] 502080.0
!K [nm^-2kjmol] 0.0
R0 [nm] 0.09572
KIND G87
&END BOND
&BOND
ATOMS OW MW
!K [nm^-2kjmol] 753120.0
K [nm^-2kjmol] 00.0
R0 [nm] 0.01750
KIND G87
&END BOND
&BOND
ATOMS HW HW
K [nm^-2kjmol] 000000.0
R0 [nm] 0.01750
KIND G87
&END BOND
&BOND
ATOMS HW MW
K [nm^-2kjmol] 000000.0
R0 [nm] 0.01750
KIND G87
&END BOND
&NONBONDED
&LENNARD-JONES
ATOMS OW OW
EPSILON [kjmol] 6.36386E-01
SIGMA [nm] 3.15061E-01
&END LENNARD-JONES
&LENNARD-JONES
ATOMS OW MW
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS OW HW
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS HW HW
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS HW MW
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&LENNARD-JONES
ATOMS MW MW
EPSILON 0.0
SIGMA 0.1
&END LENNARD-JONES
&END NONBONDED
&CHARGE
ATOM OW
! CHARGE -1.04
CHARGE 0.000000
&END CHARGE
&CHARGE
ATOM HW
CHARGE 0.52
!CHARGE 0.00
&END CHARGE
&CHARGE
ATOM MW
CHARGE -1.04
! CHARGE 0.00
&END CHARGE
&END FORCEFIELD
&POISSON
&EWALD
!EWALD_TYPE spme
EWALD_TYPE NONE
ALPHA .44
GMAX 24
O_SPLINE 6
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC 40.000 40.000 40.000
&END CELL
&TOPOLOGY
COORD_FILE_NAME ../sample_pdb/water_tip4p.pdb
COORDINATE pdb
CONN_FILE_NAME ../sample_psf/water_tip4p.psf
CONN_FILE_FORMAT PSF
&END TOPOLOGY
&KIND OW
ELEMENT O
MASS 15.99940
&END KIND
&KIND HW
ELEMENT H
MASS 1.00800
&END KIND
&KIND MW
ELEMENT H
MASS 0.00000
&END KIND
&PRINT
&ATOMIC_COORDINATES SILENT
&END ATOMIC_COORDINATES
&MOLECULES ON
&END MOLECULES
&END PRINT
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT tip4p
RUN_TYPE MD
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE
STEPS 200
TIMESTEP 1.0
TEMPERATURE 100
&END MD
&CONSTRAINT
&G3X3
INTERMOLECULAR FALSE
DISTANCES 1.89035 1.89035 3.0869
ATOMS 1 2 3
MOLECULE 1
&END G3X3
&VIRTUAL_SITE
INTERMOLECULAR FALSE
ATOMS 4 2 1 3
PARAMETERS 0.128012065 0.128012065
MOLECULE 1
&END VIRTUAL_SITE
&END CONSTRAINT
&PRINT
&FORCES ON
&EACH
MD 2
&END EACH
FORMAT XYZ
&END FORCES
&VELOCITIES ON
&EACH
MD 2
&END EACH
FORMAT XYZ
&END VELOCITIES
&END PRINT
&END MOTION

10
tests/Fist/sample_pdb/water_tip4p.pdb Normal file
View file

@ -0,0 +1,10 @@
CRYST1 0.000 0.000 0.000 90.00 90.00 90.00 P 1 1
ATOM 1 OW TIP4X 0 4.257 3.685 7.048 1.00 0.00 O
ATOM 2 HW TIP4X 0 3.505 3.065 6.819 1.00 0.00 H
ATOM 3 HW TIP4X 0 4.164 4.529 6.518 1.00 0.00 H
ATOM 4 MW TIP4X 0 4.149 3.714 6.951 1.00 0.00 H
ATOM 5 OW TIP4X 1 6.650 2.043 6.264 1.00 0.00 O
ATOM 6 HW TIP4X 1 5.872 2.583 6.586 1.00 0.00 H
ATOM 7 HW TIP4X 1 7.297 1.910 7.015 1.00 0.00 H
ATOM 8 MW TIP4X 1 6.633 2.095 6.401 1.00 0.00 H
END

41
tests/Fist/sample_psf/water_tip4p.psf Normal file
View file

@ -0,0 +1,41 @@
PSF EXT
1 !NTITLE
CP2K generated DUMP of connectivity
8 !NATOM
1 MOL1 1 TIP OW OW 0.000000 15.999 0
2 MOL1 1 TIP HW HW 0.520000 1.008 0
3 MOL1 1 TIP HW HW 0.520000 1.008 0
4 MOL1 1 TIP MW MW -1.040000 0.000 0
5 MOL1 2 TIP OW OW 0.000000 15.999 0
6 MOL1 2 TIP HW HW 0.520000 1.008 0
7 MOL1 2 TIP HW HW 0.520000 1.008 0
8 MOL1 2 TIP MW MW -1.040000 0.000 0
10 !NBOND
1 2 1 3 5 6 5 7
4 2 4 3 8 6 8 7
4 1 8 5
2 !NTHETA
2 1 3 6 5 7
0 !NPHI
0 !NIMPHI
0 !NDON
0 !NACC
0 !NNB

2
tests/SE/regtest-2/TEST_FILES_RESET
View file

@ -139,3 +139,5 @@ md.inp
scan_zn_1w_cp2k.inp
# numerical instabilities
hcl_ot.inp
# numerics
hcl_ot.inp