mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 14:15:19 -04:00
Preliminary framework for adaptive QM/MM (Noam and Teo) + 1 regtest
svn-origin-rev: 12203
This commit is contained in:
Teodoro Laino
parent
32abab9bb9
commit
70a54851ac
17 changed files with 1813 additions and 47 deletions
|
|
@ -431,6 +431,7 @@ OBJECTS_GENERIC =\
|
|||
pw_types.o\
|
||||
qmmm_elpot.o\
|
||||
qmmm_ff_fist.o\
|
||||
qmmm_force_mixing.o\
|
||||
qmmm_gaussian_data.o\
|
||||
qmmm_gaussian_init.o\
|
||||
qmmm_gaussian_input.o\
|
||||
|
|
|
|||
|
|
@ -68,7 +68,7 @@ MODULE cell_types
|
|||
! *** Public subroutines ***
|
||||
PUBLIC :: get_cell, get_cell_param, init_cell, read_cell,&
|
||||
write_cell, cell_create, cell_retain, cell_release,&
|
||||
cell_clone, compare_cells, parse_cell_line, set_cell_param
|
||||
cell_clone, cell_copy, compare_cells, parse_cell_line, set_cell_param
|
||||
|
||||
! *** Public functions ***
|
||||
PUBLIC :: plane_distance, pbc, real_to_scaled, scaled_to_real
|
||||
|
|
@ -131,6 +131,27 @@ CONTAINS
|
|||
cell_out%id_nr = last_cell_id
|
||||
END SUBROUTINE cell_clone
|
||||
|
||||
! *****************************************************************************
|
||||
SUBROUTINE cell_copy (cell_in, cell_out, error)
|
||||
TYPE(cell_type), POINTER :: cell_in, cell_out
|
||||
TYPE(cp_error_type), INTENT(INOUT) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'cell_copy', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
LOGICAL :: failure
|
||||
|
||||
failure=.FALSE.
|
||||
CPPrecondition(ASSOCIATED(cell_in),cp_failure_level,routineP,error,failure)
|
||||
CPPrecondition(ASSOCIATED(cell_out),cp_failure_level,routineP,error,failure)
|
||||
|
||||
cell_out%deth = cell_in%deth
|
||||
cell_out%perd = cell_in%perd
|
||||
cell_out%hmat = cell_in%hmat
|
||||
cell_out%h_inv = cell_in%h_inv
|
||||
cell_out%orthorhombic = cell_in%orthorhombic
|
||||
END SUBROUTINE cell_copy
|
||||
|
||||
! *****************************************************************************
|
||||
!> \brief Read cell info from a line (parsed from a file)
|
||||
!> \author Teodoro Laino [tlaino] - University of Zurich
|
||||
|
|
|
|||
|
|
@ -112,11 +112,10 @@ MODULE force_env_methods
|
|||
qmmm_added_chrg_forces,&
|
||||
qmmm_link_Imomm_coord,&
|
||||
qmmm_link_Imomm_forces
|
||||
USE qmmm_types, ONLY: fist_subsys,&
|
||||
qmmm_env_qm_retain,&
|
||||
qmmm_env_qm_type,&
|
||||
qmmm_links_type,&
|
||||
qs_subsys
|
||||
USE qmmm_types, ONLY: &
|
||||
fist_subsys, force_mixing_core_subsys, force_mixing_extended_subsys, &
|
||||
force_mixing_label_QM_extended, primary_subsys, qmmm_env_qm_retain, &
|
||||
qmmm_env_qm_type, qmmm_links_type, qs_subsys
|
||||
USE qmmm_util, ONLY: apply_qmmm_translate,&
|
||||
apply_qmmm_walls
|
||||
USE qs_energy, ONLY: qs_energies
|
||||
|
|
@ -677,7 +676,7 @@ CONTAINS
|
|||
force_env%in_use=use_qmmm
|
||||
force_env%qmmm_env => qmmm_env
|
||||
CALL qmmm_env_qm_retain(qmmm_env,error=error)
|
||||
force_env%virial => sub_force_env(1)%force_env%virial
|
||||
force_env%virial => sub_force_env(primary_subsys)%force_env%virial
|
||||
CALL virial_retain(force_env%virial,error=error)
|
||||
! Virial controlled through the external request
|
||||
CALL virial_set ( virial=force_env%virial,&
|
||||
|
|
@ -734,7 +733,7 @@ CONTAINS
|
|||
CALL eip_env_get(force_env%eip_env,subsys=force_env%subsys,error=error)
|
||||
CALL cp_subsys_retain(force_env%subsys, error=error)
|
||||
CASE(use_qmmm)
|
||||
subsys => force_env%sub_force_env(1)%force_env%subsys
|
||||
subsys => force_env%sub_force_env(primary_subsys)%force_env%subsys
|
||||
force_env%subsys => subsys
|
||||
CALL cp_subsys_retain(subsys,error=error)
|
||||
CASE(use_mixed_force)
|
||||
|
|
@ -794,7 +793,6 @@ CONTAINS
|
|||
!> \author Fawzi Mohamed
|
||||
! *****************************************************************************
|
||||
RECURSIVE SUBROUTINE qmmm_energy_and_forces(force_env,calc_force,error)
|
||||
|
||||
TYPE(force_env_type), POINTER :: force_env
|
||||
LOGICAL, INTENT(IN), OPTIONAL :: calc_force
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
|
@ -802,6 +800,90 @@ CONTAINS
|
|||
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_energy_and_forces', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: ip, isubf
|
||||
INTEGER, POINTER :: cur_indices(:), cur_labels(:)
|
||||
LOGICAL :: failure
|
||||
TYPE(cp_subsys_type), POINTER :: subsys_primary, &
|
||||
subsys_qmmm_core, &
|
||||
subsys_qmmm_extended
|
||||
TYPE(particle_type), DIMENSION(:), &
|
||||
POINTER :: particles_primary, &
|
||||
particles_qmmm_core, &
|
||||
particles_qmmm_extended
|
||||
TYPE(section_vals_type), POINTER :: force_env_section
|
||||
|
||||
! Possibly translate the system
|
||||
|
||||
CALL apply_qmmm_translate(force_env, error)
|
||||
|
||||
DO isubf=1, SIZE(force_env%sub_force_env)
|
||||
CALL qmmm_energy_and_forces_low(force_env%sub_force_env(isubf)%force_env,&
|
||||
calc_force,error)
|
||||
END DO
|
||||
|
||||
IF (SIZE(force_env%sub_force_env) == 1) THEN
|
||||
! see force_env_get for potential energy to copy
|
||||
! copy forces, energy, virial up
|
||||
CONTINUE
|
||||
ELSE IF (SIZE(force_env%sub_force_env) == 2) THEN
|
||||
! force mixing happens here
|
||||
! energy will just inherit from extended sub_force_env in force_env_get()
|
||||
|
||||
! get forces from subsys of each sub force env
|
||||
CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env,&
|
||||
subsys=subsys_qmmm_core,force_env_section=force_env_section,error=error)
|
||||
CALL force_env_get(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,&
|
||||
subsys=subsys_qmmm_extended,error=error)
|
||||
|
||||
CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%RESTART_INFO%INDICES",i_vals=cur_indices,error=error)
|
||||
CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%RESTART_INFO%LABELS",i_vals=cur_labels,error=error)
|
||||
|
||||
particles_qmmm_extended => subsys_qmmm_extended%particles%els
|
||||
particles_qmmm_core => subsys_qmmm_core%particles%els
|
||||
DO ip=1,SIZE(cur_indices)
|
||||
IF (cur_labels(ip) >= force_mixing_label_QM_extended) THEN ! this is a QM atom
|
||||
! copy (QM) force from extended calculation
|
||||
particles_qmmm_core(cur_indices(ip))%f=particles_qmmm_extended(cur_indices(ip))%f
|
||||
ENDIF
|
||||
END DO
|
||||
|
||||
IF (force_mixing_core_subsys /= primary_subsys) THEN
|
||||
CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,&
|
||||
subsys=subsys_primary,force_env_section=force_env_section,error=error)
|
||||
particles_primary => subsys_primary%particles%els
|
||||
DO ip=1,SIZE(particles_qmmm_core)
|
||||
particles_primary(ip)%f=particles_qmmm_core(ip)%f
|
||||
END DO
|
||||
ENDIF
|
||||
|
||||
ELSE
|
||||
CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,&
|
||||
routineP,"size(sub_force_env) neither 1 (conventional QM/MM) nor "//&
|
||||
"2 (force-mixing)! "//&
|
||||
CPSourceFileRef,&
|
||||
error,failure)
|
||||
ENDIF
|
||||
|
||||
END SUBROUTINE qmmm_energy_and_forces
|
||||
|
||||
! *****************************************************************************
|
||||
!> \brief calculates the qm/mm energy and forces
|
||||
!> \param calc_force if also the forces should be calculated
|
||||
!> \param error variable to control error logging, stopping,...
|
||||
!> see module cp_error_handling
|
||||
!> \par History
|
||||
!> 05.2004 created [fawzi]
|
||||
!> \author Fawzi Mohamed
|
||||
! *****************************************************************************
|
||||
RECURSIVE SUBROUTINE qmmm_energy_and_forces_low(force_env,calc_force,error)
|
||||
|
||||
TYPE(force_env_type), POINTER :: force_env
|
||||
LOGICAL, INTENT(IN), OPTIONAL :: calc_force
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_energy_and_forces_low', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: ip, output_unit
|
||||
INTEGER, DIMENSION(:), POINTER :: qm_atom_index
|
||||
LOGICAL :: calculate_forces, failure, &
|
||||
|
|
@ -849,9 +931,6 @@ CONTAINS
|
|||
END IF
|
||||
CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure)
|
||||
|
||||
! Possibly translate the system
|
||||
CALL apply_qmmm_translate(force_env, error)
|
||||
|
||||
particles_mm => subsys_mm%particles%els
|
||||
particles_qm => subsys_qm%particles%els
|
||||
|
||||
|
|
@ -925,7 +1004,7 @@ CONTAINS
|
|||
END IF
|
||||
CALL cp_print_key_finished_output(output_unit,logger,force_env_section,&
|
||||
"QMMM%PRINT%DERIVATIVES",error=error)
|
||||
END SUBROUTINE qmmm_energy_and_forces
|
||||
END SUBROUTINE qmmm_energy_and_forces_low
|
||||
|
||||
! *****************************************************************************
|
||||
!> \brief ****f* force_env_methods/mixed_energy_forces [1.0]
|
||||
|
|
|
|||
|
|
@ -60,8 +60,12 @@ MODULE force_env_types
|
|||
USE mixed_environment_types, ONLY: get_mixed_env,&
|
||||
mixed_env_release,&
|
||||
mixed_environment_type
|
||||
USE qmmm_types, ONLY: qmmm_env_qm_release,&
|
||||
qmmm_env_qm_type
|
||||
USE qmmm_types, ONLY: fist_subsys,&
|
||||
force_mixing_extended_subsys,&
|
||||
primary_subsys,&
|
||||
qmmm_env_qm_release,&
|
||||
qmmm_env_qm_type,&
|
||||
qs_subsys
|
||||
USE qs_energy_types, ONLY: qs_energy_type
|
||||
USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env
|
||||
USE qs_environment_types, ONLY: get_qs_env,&
|
||||
|
|
@ -355,6 +359,7 @@ CONTAINS
|
|||
CHARACTER(len=*), PARAMETER :: routineN = 'force_env_get', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: cur_subsys
|
||||
LOGICAL :: failure
|
||||
REAL(KIND=dp) :: eip_kinetic_energy, &
|
||||
eip_potential_energy, &
|
||||
|
|
@ -426,26 +431,38 @@ CONTAINS
|
|||
CPPrecondition(.NOT.PRESENT(kinetic_energy),cp_failure_level,routineP,error,failure)
|
||||
CASE (use_qmmm)
|
||||
IF (PRESENT(cell)) THEN
|
||||
CALL force_env_get(force_env%sub_force_env(1)%force_env,&
|
||||
CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,&
|
||||
cell=cell,&
|
||||
error=error)
|
||||
ENDIF
|
||||
IF (PRESENT(cell_ref)) THEN
|
||||
CALL force_env_get(force_env%sub_force_env(1)%force_env,&
|
||||
CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,&
|
||||
cell_ref=cell_ref,&
|
||||
error=error)
|
||||
ENDIF
|
||||
IF (PRESENT(use_ref_cell)) use_ref_cell = .FALSE.
|
||||
IF (PRESENT(kinetic_energy)) THEN
|
||||
CALL force_env_get(force_env%sub_force_env(1)%force_env,&
|
||||
CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,&
|
||||
kinetic_energy=kinetic_energy,&
|
||||
error=error)
|
||||
END IF
|
||||
IF (PRESENT(potential_energy)) THEN
|
||||
CALL force_env_get(force_env%sub_force_env(1)%force_env,&
|
||||
! get the underlying energies from primary subsys. This is the only subsys
|
||||
! for conventional QM/MM, and force-mixing knows to put relevant energy there.
|
||||
IF (SIZE(force_env%sub_force_env) == 1) THEN
|
||||
cur_subsys = primary_subsys
|
||||
ELSE IF (SIZE(force_env%sub_force_env) == 2) THEN
|
||||
cur_subsys = force_mixing_extended_subsys
|
||||
ELSE
|
||||
CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,&
|
||||
routineP,"bad number of sub_force_envs getting qmmm potential energy "//&
|
||||
CPSourceFileRef,&
|
||||
error,failure)
|
||||
ENDIF
|
||||
CALL force_env_get(force_env%sub_force_env(cur_subsys)%force_env%sub_force_env(fist_subsys)%force_env,&
|
||||
potential_energy=penergy_mm,&
|
||||
error=error)
|
||||
CALL force_env_get(force_env%sub_force_env(2)%force_env,&
|
||||
CALL force_env_get(force_env%sub_force_env(cur_subsys)%force_env%sub_force_env(qs_subsys)%force_env,&
|
||||
potential_energy=penergy_qm,&
|
||||
error=error)
|
||||
potential_energy = penergy_qm+penergy_mm
|
||||
|
|
|
|||
|
|
@ -36,7 +36,8 @@ MODULE input_cp2k_qmmm
|
|||
section_type
|
||||
USE input_val_types, ONLY: char_t,&
|
||||
integer_t,&
|
||||
real_t
|
||||
real_t, &
|
||||
logical_t
|
||||
USE kinds, ONLY: dp
|
||||
USE string_utilities, ONLY: s2a
|
||||
#include "cp_common_uses.h"
|
||||
|
|
@ -189,6 +190,11 @@ CONTAINS
|
|||
CALL section_add_keyword(section,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
! NB: remember to create these
|
||||
CALL create_qmmm_force_mixing_section(subsection, error)
|
||||
CALL section_add_subsection(section,subsection,error=error)
|
||||
CALL section_release(subsection,error=error)
|
||||
|
||||
CALL create_qmmm_qm_kinds(subsection,error)
|
||||
CALL section_add_subsection(section,subsection,error=error)
|
||||
CALL section_release(subsection,error=error)
|
||||
|
|
@ -278,6 +284,174 @@ CONTAINS
|
|||
|
||||
END SUBROUTINE create_qmmm_mm_kinds
|
||||
|
||||
! *****************************************************************************
|
||||
!> \brief Input section to create FORCE_MIXING sections
|
||||
!> \param error variable to control error logging, stopping,...
|
||||
!> see module cp_error_handling
|
||||
!> \author noam
|
||||
! *****************************************************************************
|
||||
SUBROUTINE create_qmmm_force_mixing_section(section, error)
|
||||
TYPE(section_type), POINTER :: section
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: &
|
||||
routineN = 'create_qmmm_force_mixing_section', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
LOGICAL :: failure
|
||||
TYPE(keyword_type), POINTER :: keyword
|
||||
TYPE(section_type), POINTER :: link_subsection, print_key, &
|
||||
subsection
|
||||
|
||||
failure=.FALSE.
|
||||
NULLIFY(keyword)
|
||||
CPPrecondition(.NOT.ASSOCIATED(section),cp_failure_level,routineP,error,failure)
|
||||
CALL section_create(section,name="FORCE_MIXING",&
|
||||
description="This section defines parameters for force-mixing based QM/MM,"//&
|
||||
" which actually does two conventional QM/MM calculations, on a small "//&
|
||||
" and a large QM regions, and combines the MM forces from one and QM "//&
|
||||
" forces from the other to create a complete set of forces. Energy is "//&
|
||||
" not conserved (although the energy from the large QM region is saved) "//&
|
||||
" so a proper thermostat (i.e. massive, and able to handle dissipation) must be used.",&
|
||||
n_keywords=5, n_subsections=3, repeats=.FALSE., required=.FALSE.,&
|
||||
error=error)
|
||||
|
||||
CALL keyword_create(keyword, name="_SECTION_PARAMETERS_", &
|
||||
description="Enables force-mixing", &
|
||||
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="R_CORE",&
|
||||
description="Specify the inner and outer radii of core QM region. "//&
|
||||
"All atoms within this distance (hysteretically) of any atoms listed in "//&
|
||||
"ATOM_INDEX will be core QM atoms in the force-mixing calculation.", &
|
||||
usage="R_CORE <real> <real>",n_var=2,type_of_var=real_t,&
|
||||
default_r_vals=(/ cp_unit_to_cp2k(0.0_dp,"angstrom",error=error),&
|
||||
cp_unit_to_cp2k(0.0_dp,"angstrom",error=error) /),&
|
||||
unit_str="angstrom",error=error)
|
||||
CALL section_add_keyword(section,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
CALL keyword_create(keyword, name="R_QM",&
|
||||
description="Specify the inner and outer radii of extended QM region. "//&
|
||||
"All atoms within this distance (hysteretically) of any atoms listed in "//&
|
||||
"ATOM_INDEX will be QM atoms in the force-mixing calculation.", &
|
||||
usage="R_QM <real> <real>",n_var=2,type_of_var=real_t,&
|
||||
default_r_vals=(/ cp_unit_to_cp2k(0.5_dp,"angstrom",error=error),&
|
||||
cp_unit_to_cp2k(1.0_dp,"angstrom",error=error) /),&
|
||||
unit_str="angstrom",error=error)
|
||||
CALL section_add_keyword(section,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
CALL keyword_create(keyword, name="QM_EXTENDED_SEED_IS_ONLY_CORE_LIST",&
|
||||
description="Makes the extended QM zone be defined hysterestically "//&
|
||||
"by distance from QM core list (i.e. atoms specified explicitly by "//&
|
||||
"user) instead of from full QM core region (specified by user + hysteretic "//&
|
||||
"selection + unbreakable bonds)",&
|
||||
usage="QM_EXTENDED_SEED_IS_ONLY_CORE_LIST <logical>",n_var=1,type_of_var=logical_t,&
|
||||
default_l_val=.FALSE., repeats = .FALSE., required = .FALSE., error=error)
|
||||
CALL section_add_keyword(section,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
CALL keyword_create(keyword, name="R_BUF",&
|
||||
description="Specify the inner and outer radii of buffer region. "//&
|
||||
"All atoms within this distance (hysteretically) of any QM atoms "//&
|
||||
"will be buffer atoms in the force-mixing calculation.", &
|
||||
usage="R_BUF <real> <real>",n_var=2,type_of_var=real_t,&
|
||||
default_r_vals=(/ cp_unit_to_cp2k(0.5_dp,"angstrom",error=error),&
|
||||
cp_unit_to_cp2k(1.0_dp,"angstrom",error=error) /),&
|
||||
unit_str="angstrom",error=error)
|
||||
CALL section_add_keyword(section,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
CALL keyword_create(keyword, name="QM_KIND_ELEMENT_MAPPING",&
|
||||
description="Mapping from elements to QM_KINDs for adaptivity",&
|
||||
usage="QM_KIND_ELEMENT_MAPPING {El} {QM_KIND}",&
|
||||
n_var=2,type_of_var=char_t,repeats=.TRUE.,required=.TRUE.,error=error)
|
||||
CALL section_add_keyword(section,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
CALL keyword_create(keyword, name="MAX_N_QM",&
|
||||
description="Maximum number of QM atoms, for detection of runaway adaptive selection.",&
|
||||
usage="MAX_N_QM int",default_i_val=300,&
|
||||
n_var=1,type_of_var=integer_t,repeats=.FALSE.,required=.FALSE.,error=error)
|
||||
CALL section_add_keyword(section,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
|
||||
![NB] also need a list?
|
||||
![NB] maybe not list+links , but some sort of link template
|
||||
![NB] also, breakable bonds?
|
||||
! BUFFER_LINKS subsection
|
||||
NULLIFY(subsection)
|
||||
CALL section_create(subsection,name="BUFFER_LINKS",&
|
||||
description="Information about possible links for the buffer QM/MM calculation."//&
|
||||
"Not used yet - need to implement buffer selection by atom and walking of connectivity data.",&
|
||||
n_keywords=0, n_subsections=1, repeats=.TRUE., required=.TRUE.,&
|
||||
error=error)
|
||||
|
||||
NULLIFY(link_subsection)
|
||||
CALL create_qmmm_link_section(link_subsection,error)
|
||||
CALL section_add_subsection(subsection,link_subsection,error=error)
|
||||
CALL section_release(link_subsection,error=error)
|
||||
|
||||
CALL section_add_subsection(section,subsection,error=error)
|
||||
CALL section_release(subsection,error=error)
|
||||
|
||||
! RESTART_INFO subsection
|
||||
NULLIFY(subsection)
|
||||
CALL section_create(subsection,name="RESTART_INFO",&
|
||||
description="This section provides information about old force-mixing indices and labels, "//&
|
||||
"for restarts.", &
|
||||
n_keywords=2, n_subsections=0, repeats=.FALSE., required=.TRUE.,&
|
||||
error=error)
|
||||
|
||||
CALL keyword_create(keyword, name="INDICES",&
|
||||
description="Indices of atoms in previous step QM regions.", &
|
||||
usage="INDICES 1 2 ...", &
|
||||
n_var=-1,type_of_var=integer_t,repeats=.TRUE.,required=.TRUE.,error=error)
|
||||
CALL section_add_keyword(subsection,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
CALL keyword_create(keyword, name="LABELS",&
|
||||
description="Labels of atoms in previous step QM regions.", &
|
||||
usage="LABELS 1 1 ...", &
|
||||
n_var=-1,type_of_var=integer_t,repeats=.TRUE.,required=.TRUE.,error=error)
|
||||
CALL section_add_keyword(subsection,keyword,error=error)
|
||||
CALL keyword_release(keyword,error=error)
|
||||
|
||||
CALL section_add_subsection(section,subsection,error=error)
|
||||
CALL section_release(subsection,error=error)
|
||||
|
||||
! PRINT subsection, with keys for neighbor list
|
||||
CALL section_create(subsection,name="print",&
|
||||
description="Section of possible print options in FORCE_MIXING.",&
|
||||
n_keywords=0, n_subsections=2, repeats=.FALSE., required=.TRUE.,&
|
||||
error=error)
|
||||
NULLIFY(print_key)
|
||||
CALL cp_print_key_section_create(print_key,"SUBCELL",&
|
||||
description="Activates the printing of the subcells used for the"//&
|
||||
"generation of neighbor lists.", unit_str="angstrom",&
|
||||
print_level=high_print_level,filename="__STD_OUT__",&
|
||||
error=error)
|
||||
CALL section_add_subsection(subsection,print_key,error=error)
|
||||
CALL section_release(print_key,error=error)
|
||||
|
||||
CALL cp_print_key_section_create(print_key,"NEIGHBOR_LISTS",&
|
||||
description="Activates the printing of the neighbor lists used"//&
|
||||
" for the hysteretic region calculations.", &
|
||||
print_level=high_print_level,filename="",unit_str="angstrom",&
|
||||
error=error)
|
||||
CALL section_add_subsection(subsection,print_key,error=error)
|
||||
CALL section_release(print_key,error=error)
|
||||
|
||||
CALL section_add_subsection(section,subsection,error=error)
|
||||
CALL section_release(subsection,error=error)
|
||||
|
||||
|
||||
END SUBROUTINE create_qmmm_force_mixing_section
|
||||
|
||||
! *****************************************************************************
|
||||
!> \brief Input section to create QM kinds sections
|
||||
!> \param error variable to control error logging, stopping,...
|
||||
|
|
|
|||
|
|
@ -112,7 +112,8 @@ MODULE input_cp2k_restarts
|
|||
PUBLIC :: write_restart,&
|
||||
update_subsys,&
|
||||
section_velocity_val_set,&
|
||||
write_restart_header
|
||||
write_restart_header,&
|
||||
update_input
|
||||
|
||||
CONTAINS
|
||||
|
||||
|
|
|
|||
|
|
@ -79,6 +79,7 @@ MODULE integrator
|
|||
update_particle_set
|
||||
USE physcon, ONLY: bohr,&
|
||||
femtoseconds
|
||||
USE qmmm_main, ONLY: qmmm_update_force_mixing_env
|
||||
USE qmmm_util, ONLY: apply_qmmm_walls_reflective
|
||||
USE reftraj_types, ONLY: reftraj_type
|
||||
USE reftraj_util, ONLY: compute_msd_reftraj
|
||||
|
|
@ -255,6 +256,9 @@ CONTAINS
|
|||
DEALLOCATE (pos,STAT=stat)
|
||||
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
|
||||
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
|
||||
! Update forces
|
||||
CALL force_env_calc_energy_force(force_env,error=error)
|
||||
|
||||
|
|
@ -461,6 +465,9 @@ CONTAINS
|
|||
deallocate_vel=.FALSE.
|
||||
END IF
|
||||
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
|
||||
CALL force_env_calc_energy_force(force_env,error=error)
|
||||
|
||||
! Metadynamics
|
||||
|
|
@ -627,6 +634,9 @@ CONTAINS
|
|||
CALL update_dealloc_tmp ( tmp, particle_set, shell_particle_set, &
|
||||
core_particle_set, para_env, shell_adiabatic, pos=.TRUE., error=error )
|
||||
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
|
||||
CALL force_env_calc_energy_force(force_env,error=error)
|
||||
|
||||
! Metadynamics
|
||||
|
|
@ -825,6 +835,9 @@ CONTAINS
|
|||
shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE., error=error)
|
||||
END IF
|
||||
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
|
||||
! Update forces
|
||||
CALL force_env_calc_energy_force(force_env,error=error)
|
||||
|
||||
|
|
@ -1355,6 +1368,8 @@ CONTAINS
|
|||
|
||||
! Update forces
|
||||
CALL force_env_set_cell(force_env, cell, error)
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
CALL force_env_calc_energy_force(force_env,error=error)
|
||||
|
||||
! Metadynamics
|
||||
|
|
@ -1576,6 +1591,9 @@ CPSourceFileRef,&
|
|||
CALL force_env_set_cell(force_env, cell, error)
|
||||
END IF
|
||||
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
|
||||
! Task to perform on the reference trajectory
|
||||
! Compute energy and forces
|
||||
CALL force_env_calc_energy_force(force_env,eval_energy_forces=reftraj_env%info%eval_EF,error=error)
|
||||
|
|
@ -1797,6 +1815,8 @@ CPSourceFileRef,&
|
|||
|
||||
! Update forces (and stress)
|
||||
CALL force_env_set_cell(force_env, cell, error)
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
CALL force_env_calc_energy_force(force_env, error=error)
|
||||
|
||||
! Metadynamics
|
||||
|
|
@ -2059,6 +2079,8 @@ CPSourceFileRef,&
|
|||
|
||||
! Update forces
|
||||
CALL force_env_set_cell(force_env, cell, error)
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
CALL force_env_calc_energy_force(force_env, error=error)
|
||||
|
||||
! Metadynamics
|
||||
|
|
@ -2357,6 +2379,8 @@ CPSourceFileRef,&
|
|||
|
||||
! Update forces
|
||||
CALL force_env_set_cell(force_env, cell, error)
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
CALL force_env_calc_energy_force(force_env,error=error)
|
||||
|
||||
! Metadynamics
|
||||
|
|
@ -2579,6 +2603,9 @@ CPSourceFileRef,&
|
|||
particle_set_respa(iparticle)%r = particle_set(iparticle)%r
|
||||
END DO
|
||||
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
! Not checked, and probably won't work, with RESPA.
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
! Update forces
|
||||
CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env,error=error)
|
||||
|
||||
|
|
@ -2611,6 +2638,9 @@ CPSourceFileRef,&
|
|||
|
||||
! Multiple time step (second part)
|
||||
! Compute forces for respa force_env
|
||||
![ADAPT] update input structure with new coordinates, make new labels
|
||||
! Not checked, and probably won't work, with RESPA.
|
||||
CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error)
|
||||
CALL force_env_calc_energy_force(force_env,error=error)
|
||||
|
||||
! Metadynamics
|
||||
|
|
|
|||
735
src/qmmm_force_mixing.F
Normal file
735
src/qmmm_force_mixing.F
Normal file
|
|
@ -0,0 +1,735 @@
|
|||
! *****************************************************************************
|
||||
!> \brief Routines used for force-mixing QM/MM calculations
|
||||
!> \par History
|
||||
!> 2.2012 created [noam]
|
||||
!> \author Noam Bernstein
|
||||
! *****************************************************************************
|
||||
MODULE qmmm_force_mixing
|
||||
USE atomic_kind_types, ONLY: atomic_kind_type,&
|
||||
get_atomic_kind,&
|
||||
set_atomic_kind
|
||||
USE cell_types, ONLY: cell_type,&
|
||||
pbc
|
||||
USE cp_subsys_types, ONLY: cp_subsys_get,&
|
||||
cp_subsys_type
|
||||
USE f77_blas
|
||||
USE fist_neighbor_list_types, ONLY: fist_neighbor_deallocate,&
|
||||
fist_neighbor_type
|
||||
USE fist_neighbor_lists, ONLY: build_fist_neighbor_lists
|
||||
USE force_env_types, ONLY: force_env_get,&
|
||||
force_env_type
|
||||
USE input_section_types, ONLY: &
|
||||
section_vals_add_values, section_vals_duplicate, section_vals_get, &
|
||||
section_vals_get_subs_vals, section_vals_get_subs_vals3, &
|
||||
section_vals_remove_values, section_vals_type, section_vals_val_get, &
|
||||
section_vals_val_set
|
||||
USE kinds, ONLY: default_string_length,&
|
||||
dp
|
||||
USE memory_utilities, ONLY: reallocate
|
||||
USE mol_new_list_types, ONLY: mol_new_list_type
|
||||
USE molecule_types_new, ONLY: molecule_type
|
||||
USE particle_list_types, ONLY: particle_list_type
|
||||
USE particle_types, ONLY: particle_type
|
||||
USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
|
||||
USE qmmm_types, ONLY: force_mixing_label_QM_core,&
|
||||
force_mixing_label_QM_core_list,&
|
||||
force_mixing_label_QM_extended,&
|
||||
force_mixing_label_buffer,&
|
||||
force_mixing_label_none,&
|
||||
force_mixing_label_termination
|
||||
#include "cp_common_uses.h"
|
||||
|
||||
IMPLICIT NONE
|
||||
PRIVATE
|
||||
|
||||
LOGICAL, PRIVATE :: debug_this_module=.FALSE.
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_force_mixing'
|
||||
|
||||
PUBLIC :: setup_force_mixing_qmmm_sections, update_force_mixing_labels
|
||||
|
||||
|
||||
!***
|
||||
CONTAINS
|
||||
|
||||
! *****************************************************************************
|
||||
!> \param error variable to control error logging, stopping,...
|
||||
!> see module cp_error_handling
|
||||
!> \par History
|
||||
!> 02.2012 created [noam]
|
||||
!> \author Noam Bernstein
|
||||
! *****************************************************************************
|
||||
SUBROUTINE update_force_mixing_labels(force_env, subsys, qmmm_section, labels_changed, error)
|
||||
TYPE(force_env_type), POINTER :: force_env
|
||||
TYPE(cp_subsys_type), POINTER :: subsys
|
||||
TYPE(section_vals_type), POINTER :: qmmm_section
|
||||
LOGICAL, OPTIONAL :: labels_changed
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'update_force_mixing_labels', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: i_rep_section, i_rep_val, ip, im, &
|
||||
n_new, n_rep_section, &
|
||||
n_rep_val, natoms
|
||||
INTEGER, ALLOCATABLE :: new_full_labels(:), &
|
||||
orig_full_labels(:)
|
||||
INTEGER, POINTER :: cur_indices(:), &
|
||||
cur_labels(:), &
|
||||
mm_index_entry(:), &
|
||||
new_indices(:), new_labels(:)
|
||||
LOGICAL :: explicit
|
||||
REAL(dp), ALLOCATABLE :: nearest_dist(:)
|
||||
REAL(dp), POINTER :: r_buf(:), r_core(:), r_qm(:)
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
TYPE(fist_neighbor_type), POINTER :: nlist
|
||||
TYPE(mol_new_list_type), POINTER :: molecules
|
||||
TYPE(molecule_type), DIMENSION(:), &
|
||||
POINTER :: molecule_set
|
||||
TYPE(particle_list_type), POINTER :: particles
|
||||
TYPE(particle_type), DIMENSION(:), &
|
||||
POINTER :: particle_set
|
||||
TYPE(section_vals_type), POINTER :: force_mixing_section, &
|
||||
qm_kind_section, &
|
||||
restart_section
|
||||
integer :: QM_extended_seed_min_label_val
|
||||
logical :: QM_extended_seed_is_core_list
|
||||
integer, pointer :: breakable_bonds(:,:), broken_bonds(:)
|
||||
integer, pointer :: molecule_of_atom(:)
|
||||
integer :: max_n_qm
|
||||
|
||||
if (debug_this_module) print *, "BOB starting update_force_mixing_labels"
|
||||
! get cur indices, labels
|
||||
force_mixing_section => section_vals_get_subs_vals3(qmmm_section,"FORCE_MIXING",error=error)
|
||||
CALL get_force_mixing_indices(force_mixing_section, cur_indices, cur_labels, error=error)
|
||||
if (debug_this_module) print *, "BOB got cur_indices ",cur_indices(1:n_new)
|
||||
if (debug_this_module) print *, "BOB got cur_labels ",cur_labels(1:n_new)
|
||||
|
||||
! read from input
|
||||
![NB] breakable bonds will come from here, too
|
||||
NULLIFY(r_core, r_qm, r_buf)
|
||||
CALL section_vals_val_get(force_mixing_section,"R_CORE",r_vals=r_core,error=error)
|
||||
CALL section_vals_val_get(force_mixing_section,"R_QM",r_vals=r_qm,error=error)
|
||||
CALL section_vals_val_get(force_mixing_section,"QM_EXTENDED_SEED_IS_ONLY_CORE_LIST",&
|
||||
l_val=QM_extended_seed_is_core_list,error=error)
|
||||
CALL section_vals_val_get(force_mixing_section,"R_BUF",r_vals=r_buf,error=error)
|
||||
CALL section_vals_val_get(force_mixing_section,"MAX_N_QM",i_val=max_n_qm,error=error)
|
||||
![NB] need to read real list from input
|
||||
! should be 2xN_bb integer arrays, with (1,:) indices of inside atoms, and (2,:) indices of outside atoms
|
||||
! maybe also breakable_bond_types, with _atomic numbers_ of inside/outside atoms?
|
||||
! separate lists for core/buffer?
|
||||
allocate(breakable_bonds(2,1))
|
||||
breakable_bonds(1:2,1) = 0
|
||||
|
||||
! get particles, molecules
|
||||
NULLIFY(particles,molecules)
|
||||
CALL cp_subsys_get(subsys=subsys, particles=particles,molecules_new=molecules, error=error)
|
||||
particle_set => particles%els
|
||||
molecule_set => molecules%els
|
||||
|
||||
natoms = SIZE(particle_set)
|
||||
|
||||
! initialize new indices, labels, and new_full_labels
|
||||
NULLIFY(new_indices, new_labels)
|
||||
CALL reallocate(new_indices,1,SIZE(cur_indices))
|
||||
CALL reallocate(new_labels,1,SIZE(cur_labels))
|
||||
new_indices = 0
|
||||
new_labels = force_mixing_label_none
|
||||
ALLOCATE(new_full_labels(natoms))
|
||||
new_full_labels = force_mixing_label_none
|
||||
|
||||
! neighbor list for various hysteretic distance calls
|
||||
NULLIFY(cell)
|
||||
CALL force_env_get(force_env, cell=cell, error=error)
|
||||
NULLIFY(nlist)
|
||||
CALL make_neighbor_list(force_mixing_section, subsys, cell, MAX(r_core(2), r_qm(2), r_buf(2)), nlist, error)
|
||||
|
||||
! reset core_list from QM_KIND
|
||||
NULLIFY(mm_index_entry)
|
||||
qm_kind_section => section_vals_get_subs_vals3(qmmm_section,"QM_KIND",error=error)
|
||||
CALL section_vals_get(qm_kind_section,n_repetition=n_rep_section,error=error)
|
||||
n_new = 0
|
||||
DO i_rep_section=1,n_rep_section
|
||||
CALL section_vals_val_get(qm_kind_section,"MM_INDEX",i_rep_section=i_rep_section,n_rep_val=n_rep_val,error=error)
|
||||
DO i_rep_val=1,n_rep_val
|
||||
CALL section_vals_val_get(qm_kind_section,"MM_INDEX",i_rep_section=i_rep_section,i_rep_val=i_rep_val, &
|
||||
i_vals=mm_index_entry, error=error)
|
||||
DO ip=1, SIZE(mm_index_entry)
|
||||
call add_new_label(mm_index_entry(ip), force_mixing_label_QM_core_list, n_new, new_indices, new_labels, &
|
||||
new_full_labels, max_n_qm)
|
||||
END DO ! ip
|
||||
END DO ! i_rep_val
|
||||
END DO ! i_rep_section
|
||||
|
||||
if (debug_this_module) print *, "BOB core_list new_indices ",new_indices(1:n_new)
|
||||
if (debug_this_module) print *, "BOB core_list new_labels ",new_labels(1:n_new)
|
||||
|
||||
! allocate and initialize full atom set labels for hysteretic loops
|
||||
ALLOCATE(nearest_dist(natoms))
|
||||
|
||||
! orig_full_labels is full array (natoms) with orig labels
|
||||
ALLOCATE(orig_full_labels(natoms))
|
||||
orig_full_labels = force_mixing_label_none
|
||||
orig_full_labels(cur_indices(:)) = cur_labels(:)
|
||||
|
||||
! hysteretically set QM core from QM_core_list and radii, whole molecule
|
||||
![NB] need to replace all the whole molecule stuff with pad to breakable bonds. not quite done
|
||||
! (need intra molecule bond info, which isn't available for QM molecules yet)
|
||||
|
||||
! find molecule of atom
|
||||
allocate(molecule_of_atom(natoms))
|
||||
do im=1, size(molecule_set)
|
||||
molecule_of_atom(molecule_set(im)%first_atom:molecule_set(im)%last_atom) = im
|
||||
end do
|
||||
|
||||
! add core using hysteretic selection(core_list, r_core) + unbreakable bonds
|
||||
CALL add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, &
|
||||
orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
|
||||
force_mixing_label_QM_core_list, force_mixing_label_QM_core_list, force_mixing_label_QM_core, r_core, &
|
||||
max_n_qm, .TRUE., breakable_bonds, molecule_of_atom, molecule_set, broken_bonds)
|
||||
![NB] should actually pass this back for making link sections?
|
||||
deallocate(broken_bonds)
|
||||
|
||||
if (debug_this_module) print *, "BOB core new_indices ",new_indices(1:n_new)
|
||||
if (debug_this_module) print *, "BOB core new_labels ",new_labels(1:n_new)
|
||||
|
||||
![NB] need more sophisticated QM extended, buffer rules
|
||||
|
||||
! add QM using hysteretic selection (core_list, r_qm) + unbreakable bonds
|
||||
if (debug_this_module) print *, "BOB QM_extended_seed_is_core_list ",QM_extended_seed_is_core_list
|
||||
if (QM_extended_seed_is_core_list) then
|
||||
QM_extended_seed_min_label_val = force_mixing_label_QM_core_list
|
||||
else ! QM region seed is all of core, not just core list + unbreakable bonds
|
||||
QM_extended_seed_min_label_val = force_mixing_label_QM_core
|
||||
endif
|
||||
CALL add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, &
|
||||
orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
|
||||
QM_extended_seed_min_label_val, force_mixing_label_QM_core_list, &
|
||||
force_mixing_label_QM_extended, r_qm, &
|
||||
max_n_qm, .TRUE., breakable_bonds, molecule_of_atom, molecule_set)
|
||||
|
||||
if (debug_this_module) print *, "BOB extended new_indices ",new_indices(1:n_new)
|
||||
if (debug_this_module) print *, "BOB extended new_labels ",new_labels(1:n_new)
|
||||
|
||||
! add buffer using hysteretic selection (>= QM extended, r_buf) + unbreakable bonds
|
||||
CALL add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, &
|
||||
orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
|
||||
force_mixing_label_QM_extended, force_mixing_label_QM_core_list, force_mixing_label_buffer, r_buf, &
|
||||
max_n_qm, .TRUE., breakable_bonds, molecule_of_atom, molecule_set, broken_bonds)
|
||||
![NB] should actually pass this back for making link sections?
|
||||
deallocate(broken_bonds)
|
||||
|
||||
if (debug_this_module) print *, "BOB buffer new_indices ",new_indices(1:n_new)
|
||||
if (debug_this_module) print *, "BOB buffer new_labels ",new_labels(1:n_new)
|
||||
|
||||
deallocate(molecule_of_atom, breakable_bonds)
|
||||
DEALLOCATE(nearest_dist)
|
||||
|
||||
IF (PRESENT(labels_changed)) labels_changed = ANY (new_full_labels /= orig_full_labels)
|
||||
|
||||
DEALLOCATE(orig_full_labels)
|
||||
DEALLOCATE(new_full_labels)
|
||||
|
||||
! reduce new indices, labels to actually used size
|
||||
CALL reallocate(new_indices,1,n_new)
|
||||
CALL reallocate(new_labels,1,n_new)
|
||||
|
||||
! save info in input structure
|
||||
restart_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING%RESTART_INFO", error=error)
|
||||
CALL section_vals_get(restart_section,explicit=explicit,error=error)
|
||||
IF (explicit) CALL section_vals_remove_values(restart_section,error)
|
||||
CALL section_vals_val_set(restart_section,"INDICES",i_vals_ptr=new_indices,error=error)
|
||||
CALL section_vals_val_set(restart_section,"LABELS",i_vals_ptr=new_labels,error=error)
|
||||
|
||||
DEALLOCATE(cur_indices, cur_labels)
|
||||
CALL fist_neighbor_deallocate(nlist,error)
|
||||
|
||||
![NB] perhap be controlled by some &PRINT section?
|
||||
if (force_env%para_env%ionode) then
|
||||
print *, "QMMM FORCE MIXING final count (not including links): "//&
|
||||
" N_QM core ", count(new_labels >= force_mixing_label_QM_core), &
|
||||
" N_QM buffered ", count(new_labels >= force_mixing_label_buffer)
|
||||
endif
|
||||
|
||||
END SUBROUTINE update_force_mixing_labels
|
||||
|
||||
SUBROUTINE add_new_label(ip, label, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
|
||||
integer :: ip, label, n_new
|
||||
integer, pointer :: new_indices(:), new_labels(:)
|
||||
integer :: new_full_labels(:)
|
||||
integer :: max_n_qm
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'update_force_mixing_labels', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
n_new = n_new + 1
|
||||
CALL cp_assert(n_new <= max_n_qm, cp_failure_level,cp_assertion_failed,routineP,&
|
||||
"add_new_label tried to add more atoms than allowed by &FORCE_MIXING&MAX_N_QM!"//&
|
||||
CPSourceFileRef)
|
||||
IF (n_new > SIZE(new_indices)) CALL reallocate(new_indices,1,n_new+9)
|
||||
IF (n_new > SIZE(new_labels)) CALL reallocate(new_labels,1,n_new+9)
|
||||
new_indices(n_new) = ip
|
||||
new_labels(n_new) = label
|
||||
new_full_labels(ip) = label
|
||||
END SUBROUTINE add_new_label
|
||||
|
||||
|
||||
SUBROUTINE add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, &
|
||||
orig_full_labels, new_full_labels, n_new, new_indices, new_labels, &
|
||||
seed_min_label_val, seed_max_label_val, set_label_val, r_inout, max_n_qm, &
|
||||
pad_to_breakable_bonds, breakable_bonds, molecule_of_atom, molecule_set, &
|
||||
broken_bonds)
|
||||
TYPE(fist_neighbor_type), POINTER :: nlist
|
||||
TYPE(particle_type), DIMENSION(:), &
|
||||
POINTER :: particle_set
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
REAL(dp) :: nearest_dist(:)
|
||||
INTEGER :: orig_full_labels(:), new_full_labels(:)
|
||||
integer :: n_new
|
||||
INTEGER, POINTER :: new_indices(:), new_labels(:)
|
||||
INTEGER :: seed_min_label_val, &
|
||||
seed_max_label_val, &
|
||||
set_label_val
|
||||
REAL(dp) :: r_inout(2)
|
||||
integer :: max_n_qm
|
||||
LOGICAL :: pad_to_breakable_bonds
|
||||
integer, pointer, optional :: molecule_of_atom(:), breakable_bonds(:,:), broken_bonds(:)
|
||||
integer :: bond_atom_a, bond_atom_b, bond_atom_inside, bond_atom_outside
|
||||
TYPE(molecule_type), DIMENSION(:), &
|
||||
OPTIONAL, POINTER :: molecule_set
|
||||
|
||||
integer :: natoms
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'add_layer_hysteretically', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: i_ind, im, ip, ib, ibb, ipair, &
|
||||
ipairkind, j_ind
|
||||
LOGICAL :: molec_in_inner, molec_in_outer
|
||||
LOGICAL :: i_in_new_seed, j_in_new_seed
|
||||
LOGICAL :: i_outside_new_seed, j_outside_new_seed
|
||||
LOGICAL :: breakable
|
||||
REAL(dp) :: r_ij(3), r_ij_mag
|
||||
integer :: orig_n_new
|
||||
logical :: any_additions
|
||||
integer :: n_broken_bonds
|
||||
|
||||
if (debug_this_module) print *,"BOB adding hysteretically seed ",seed_min_label_val,seed_max_label_val,&
|
||||
" set ",set_label_val," r ",r_inout
|
||||
! calculate nearest dist from each atom outside of new seed to nearest atom inside of new seed
|
||||
nearest_dist=HUGE(1.0_dp)
|
||||
! loop over pairs of all kinds in random order
|
||||
DO ipairkind=1,SIZE(nlist%neighbor_kind_pairs)
|
||||
DO ipair=1,nlist%neighbor_kind_pairs(ipairkind)%npairs
|
||||
|
||||
i_ind = nlist%neighbor_kind_pairs(ipairkind)%list(1,ipair)
|
||||
j_ind = nlist%neighbor_kind_pairs(ipairkind)%list(2,ipair)
|
||||
|
||||
i_in_new_seed = (new_full_labels(i_ind) >= seed_min_label_val .AND. new_full_labels(i_ind) <= seed_max_label_val)
|
||||
i_outside_new_seed = (new_full_labels(i_ind) < seed_min_label_val)
|
||||
j_in_new_seed = (new_full_labels(j_ind) >= seed_min_label_val .AND. new_full_labels(j_ind) <= seed_max_label_val)
|
||||
j_outside_new_seed = (new_full_labels(j_ind) < seed_min_label_val)
|
||||
|
||||
IF ((i_in_new_seed .and. j_outside_new_seed) .or. (j_in_new_seed .and. i_outside_new_seed)) THEN
|
||||
r_ij = pbc(particle_set(i_ind)%r - particle_set(j_ind)%r, cell)
|
||||
r_ij_mag = SQRT(SUM(r_ij**2))
|
||||
IF (i_in_new_seed .and. j_outside_new_seed .and. (r_ij_mag < nearest_dist(j_ind))) then
|
||||
nearest_dist(j_ind) = r_ij_mag
|
||||
ENDIF
|
||||
IF (j_in_new_seed .and. i_outside_new_seed .and. (r_ij_mag < nearest_dist(i_ind))) then
|
||||
nearest_dist(i_ind) = r_ij_mag
|
||||
ENDIF
|
||||
ENDIF
|
||||
|
||||
END DO
|
||||
END DO
|
||||
if (debug_this_module) print *,"BOB nearest_dist ", nearest_dist
|
||||
|
||||
![NB] this is whole molecule. Should be replaced with labeling of individual atoms +
|
||||
! pad_to_breakable_bonds (below), but QM molecule bond information isn't available yet
|
||||
DO im=1, SIZE(molecule_set)
|
||||
if (debug_this_module) print *, "BOB check molecule ",im," atoms ",molecule_set(im)%first_atom,molecule_set(im)%last_atom
|
||||
molec_in_inner = any(nearest_dist(molecule_set(im)%first_atom:molecule_set(im)%last_atom) <= r_inout(1))
|
||||
molec_in_outer = any(nearest_dist(molecule_set(im)%first_atom:molecule_set(im)%last_atom) <= r_inout(2))
|
||||
if (debug_this_module) print *, "BOB got molec_in_inner ", molec_in_inner, " outer ", molec_in_outer
|
||||
if (molec_in_inner) then
|
||||
do ip=molecule_set(im)%first_atom, molecule_set(im)%last_atom
|
||||
! labels are being rebuild from scratch, so never overwrite new label that's higher level
|
||||
if (new_full_labels(ip) < set_label_val) &
|
||||
call add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
|
||||
end do
|
||||
else if (molec_in_outer) then
|
||||
if (any(orig_full_labels(molecule_set(im)%first_atom:molecule_set(im)%last_atom) >= set_label_val)) then
|
||||
do ip=molecule_set(im)%first_atom, molecule_set(im)%last_atom
|
||||
! labels are being rebuild from scratch, so never overwrite new label that's higher level
|
||||
if (new_full_labels(ip) < set_label_val) &
|
||||
call add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
|
||||
end do
|
||||
endif
|
||||
endif
|
||||
end do
|
||||
if (present(broken_bonds)) CALL reallocate(broken_bonds,1,0)
|
||||
|
||||
!!![NB] enable this (and remove previous bit) once intra-QM-molecule bond information
|
||||
!!! is available
|
||||
!!! ! label new layer
|
||||
!!! ! loop over every atom
|
||||
!!! DO ip=1, SIZE(new_full_labels)
|
||||
!!! ! labels are being rebuild from scratch, so never overwrite new label that's higher level
|
||||
!!! if (new_full_labels(ip) > set_label_val) cycle
|
||||
!!!
|
||||
!!!if (debug_this_module) print *,"BOB atom ",ip," new_full_label ",new_full_labels(ip), &
|
||||
!!!" is not in seed, check for for adding"
|
||||
!!! ![NB] is this logic right? maybe it is
|
||||
!!! IF ( nearest_dist(ip) <= r_inout(1) .OR. &
|
||||
!!! (nearest_dist(ip) <= r_inout(2) .AND. orig_full_labels(ip) >= set_label_val) ) THEN
|
||||
!!! ! atom in inner radius
|
||||
!!! ! or
|
||||
!!! ! atom in outer radius and was previously set to set_label_val or higher
|
||||
!!!if (debug_this_module) print *,"BOB add label for atom ", ip
|
||||
!!! call add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
|
||||
!!! ENDIF
|
||||
!!! END DO ! ip
|
||||
!!!
|
||||
!!! if (present(broken_bonds)) n_broken_bonds = 0
|
||||
!!! if (pad_to_breakable_bonds) then
|
||||
!!!if (debug_this_module) print *,"BOB pad to breakable bonds"
|
||||
!!! CALL cp_assert(PRESENT(molecule_of_atom),cp_failure_level,cp_assertion_failed,routineP,&
|
||||
!!! "pad_to_breakable_bonds is true, but molecule_of_atom not present!"//&
|
||||
!!!CPSourceFileRef)
|
||||
!!! CALL cp_assert(PRESENT(molecule_set),cp_failure_level,cp_assertion_failed,routineP,&
|
||||
!!! "pad_to_breakable_bonds is true, but molecule_set not present!"//&
|
||||
!!!CPSourceFileRef)
|
||||
!!! CALL cp_assert(PRESENT(breakable_bonds),cp_failure_level,cp_assertion_failed,routineP,&
|
||||
!!! "whole_molecule is true, but breakable_bonds not present!"//&
|
||||
!!!CPSourceFileRef)
|
||||
!!!
|
||||
!!! orig_n_new = n_new
|
||||
!!! ! loop over all recently marked atoms
|
||||
!!! do ip=1, orig_n_new
|
||||
!!! if (new_labels(ip) /= set_label_val) cycle
|
||||
!!!
|
||||
!!!if (debug_this_module) then
|
||||
!!! print *, "BOB molecule set"
|
||||
!!! do im=1, size(molecule_set)
|
||||
!!! print *, "BOB molecule im ",im, " #bonds ",size(molecule_set(im)%molecule_kind%bond_list), " name ", &
|
||||
!!! trim(molecule_set(im)%molecule_kind%name), molecule_set(im)%molecule_kind%kind_number
|
||||
!!! do ib=1, size(molecule_set(im)%molecule_kind%bond_list)
|
||||
!!! print *, "BOB",molecule_set(im)%molecule_kind%bond_list(ib)%a + molecule_set(im)%first_atom - 1, &
|
||||
!!! molecule_set(im)%molecule_kind%bond_list(ib)%b + molecule_set(im)%first_atom - 1
|
||||
!!! end do
|
||||
!!! end do
|
||||
!!!endif
|
||||
!!!
|
||||
!!!if (debug_this_module) print *,"BOB check atom local index ",ip," full index ", new_indices(ip)
|
||||
!!! ! atom has just been set to set_label_val, check the rest of its molecule
|
||||
!!! im = molecule_of_atom(new_indices(ip))
|
||||
!!!
|
||||
!!!if (debug_this_module) print *,"in molecule ",im
|
||||
!!! any_additions = .true.
|
||||
!!! do while (any_additions)
|
||||
!!!if (debug_this_module) print *,"start loop checking for additions"
|
||||
!!! any_additions = .false.
|
||||
!!! ! loop over bonds in this molecule in random order
|
||||
!!! do ib=1, size(molecule_set(im)%molecule_kind%bond_list)
|
||||
!!! bond_atom_a = molecule_set(im)%molecule_kind%bond_list(ib)%a + molecule_set(im)%first_atom - 1
|
||||
!!! bond_atom_b = molecule_set(im)%molecule_kind%bond_list(ib)%b + molecule_set(im)%first_atom - 1
|
||||
!!!if (debug_this_module) print *,"bond ", ib, " between a,b ", bond_atom_a,bond_atom_b
|
||||
!!!
|
||||
!!! if (new_full_labels(bond_atom_a) == set_label_val .and. new_full_labels(bond_atom_b) < set_label_val) then
|
||||
!!! ! a is inside, b is outside
|
||||
!!! bond_atom_inside = bond_atom_a
|
||||
!!! bond_atom_outside = bond_atom_b
|
||||
!!!if (debug_this_module) print *," a is inside, b is outside"
|
||||
!!! else if (new_full_labels(bond_atom_b) == set_label_val .and. new_full_labels(bond_atom_a) < set_label_val) then
|
||||
!!! ! b is inside, a is outside
|
||||
!!! bond_atom_inside = bond_atom_b
|
||||
!!! bond_atom_outside = bond_atom_a
|
||||
!!!if (debug_this_module) print *," b is inside, a is outside"
|
||||
!!! else
|
||||
!!! ! inside-inside or outside-outside
|
||||
!!! cycle
|
||||
!!! endif
|
||||
!!!
|
||||
!!! ! check for breakable bonds
|
||||
!!! breakable = .false.
|
||||
!!! do ibb=1,size(breakable_bonds,2)
|
||||
!!! if ((breakable_bonds(1,ibb) == bond_atom_inside) .and. breakable_bonds(2,ibb) == bond_atom_outside) then
|
||||
!!! ! this bond is OK to break.
|
||||
!!! breakable = .true.
|
||||
!!! if (set_label_val == force_mixing_label_buffer) then ! set termination mark on outside atom b
|
||||
!!! call add_new_label(bond_atom_outside, force_mixing_label_termination, n_new, new_indices, new_labels, &
|
||||
!!! new_full_labels, max_n_qm)
|
||||
!!! endif
|
||||
!!! if (present(broken_bonds)) then ! add a new one
|
||||
!!! n_broken_bonds = n_broken_bonds + 1
|
||||
!!! CALL reallocate(broken_bonds,1,n_broken_bonds)
|
||||
!!! broken_bonds(n_broken_bonds) = ibb
|
||||
!!! endif
|
||||
!!! endif
|
||||
!!! end do ! ibb
|
||||
!!!if (debug_this_module) print *," breakable ", breakable
|
||||
!!!
|
||||
!!! if (.not. breakable) then ! add correct atom
|
||||
!!! call add_new_label(bond_atom_outside, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm)
|
||||
!!! any_additions = .true.
|
||||
!!!if (debug_this_module) print *, " not breakable, adding outside atom ", bond_atom_outside
|
||||
!!! endif
|
||||
!!!
|
||||
!!! end do ! ib
|
||||
!!! end do ! any_additions
|
||||
!!!
|
||||
!!! end do ! ip
|
||||
!!!
|
||||
!!! if (present(broken_bonds)) CALL reallocate(broken_bonds,1,n_broken_bonds)
|
||||
!!! endif ! pad_to_breakable_bonds
|
||||
|
||||
END SUBROUTINE add_layer_hysteretically
|
||||
|
||||
SUBROUTINE make_neighbor_list(force_mixing_section, subsys, cell, r_max, nlist, error)
|
||||
TYPE(section_vals_type), POINTER :: force_mixing_section
|
||||
TYPE(cp_subsys_type), POINTER :: subsys
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
REAL(dp) :: r_max
|
||||
TYPE(fist_neighbor_type), POINTER :: nlist
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'make_neighbor_list', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
CHARACTER(LEN=default_string_length) :: kind_name
|
||||
CHARACTER(LEN=default_string_length), &
|
||||
POINTER :: kind_name_a(:)
|
||||
INTEGER :: ik
|
||||
LOGICAL :: skip_kind
|
||||
REAL(dp), ALLOCATABLE :: r_max_a(:,:), r_minsq_a(:,:)
|
||||
TYPE(atomic_kind_type), POINTER :: atomic_kind
|
||||
|
||||
ALLOCATE(r_max_a(SIZE(subsys%atomic_kinds%els),SIZE(subsys%atomic_kinds%els)))
|
||||
ALLOCATE(r_minsq_a(SIZE(subsys%atomic_kinds%els),SIZE(subsys%atomic_kinds%els)))
|
||||
r_max_a = r_max
|
||||
r_minsq_a = EPSILON(1.0_dp)
|
||||
|
||||
! save kind names
|
||||
ALLOCATE(kind_name_a(SIZE(subsys%atomic_kinds%els)))
|
||||
DO ik=1, SIZE(subsys%atomic_kinds%els)
|
||||
atomic_kind => subsys%atomic_kinds%els(ik)
|
||||
CALL get_atomic_kind(atomic_kind=atomic_kind,name=kind_name)
|
||||
kind_name_a(ik) = kind_name
|
||||
END DO
|
||||
|
||||
! overwrite kind names so that none are QM, and so excluding QM-QM interactions
|
||||
! (which is not what we want) will not happen
|
||||
DO ik=1, SIZE(subsys%atomic_kinds%els)
|
||||
atomic_kind => subsys%atomic_kinds%els(ik)
|
||||
CALL get_atomic_kind(atomic_kind=atomic_kind,name=kind_name)
|
||||
! when atom is QM atom, kind_name is replaced with original
|
||||
! mm kind name, and return status is logical .TRUE.
|
||||
skip_kind=qmmm_ff_precond_only_qm(kind_name)
|
||||
CALL set_atomic_kind(atomic_kind=atomic_kind,name=kind_name)
|
||||
END DO
|
||||
|
||||
NULLIFY(nlist)
|
||||
CALL build_fist_neighbor_lists(subsys%atomic_kinds%els, subsys%particles%els, &
|
||||
cell=cell, r_max=r_max_a, r_minsq=r_minsq_a, &
|
||||
ei_scale14=1.0_dp, vdw_scale14=1.0_dp, nonbonded=nlist, &
|
||||
para_env=subsys%para_env, build_from_scratch=.TRUE., geo_check=.FALSE., &
|
||||
mm_section=force_mixing_section, error=error)
|
||||
|
||||
DEALLOCATE(r_max_a, r_minsq_a)
|
||||
|
||||
! restore kind names
|
||||
DO ik=1, SIZE(subsys%atomic_kinds%els)
|
||||
CALL set_atomic_kind(atomic_kind=atomic_kind,name=kind_name_a(ik))
|
||||
END DO
|
||||
DEALLOCATE(kind_name_a)
|
||||
|
||||
END SUBROUTINE make_neighbor_list
|
||||
|
||||
! *****************************************************************************
|
||||
!> \param error variable to control error logging, stopping,...
|
||||
!> see module cp_error_handling
|
||||
!> \par History
|
||||
!> 02.2012 created [noam]
|
||||
!> \author Noam Bernstein
|
||||
! *****************************************************************************
|
||||
SUBROUTINE setup_force_mixing_qmmm_sections(subsys, qmmm_section, qmmm_core_section, qmmm_extended_section, error)
|
||||
TYPE(cp_subsys_type), POINTER :: subsys
|
||||
TYPE(section_vals_type), POINTER :: qmmm_section, &
|
||||
qmmm_core_section, &
|
||||
qmmm_extended_section
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: &
|
||||
routineN = 'setup_force_mixing_qmmm_sections', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
CHARACTER(len=DEFAULT_STRING_LENGTH), &
|
||||
POINTER :: elem_mapping(:,:), &
|
||||
elem_mapping_entry(:)
|
||||
INTEGER :: i_rep_section_core, i_rep_section_extended, i_rep_val_core, &
|
||||
i_rep_val_extended, ielem, ip, n_elements
|
||||
INTEGER, POINTER :: cur_indices(:), cur_labels(:)
|
||||
LOGICAL :: new_element_core, &
|
||||
new_element_extended
|
||||
TYPE(particle_type), DIMENSION(:), &
|
||||
POINTER :: particles
|
||||
TYPE(section_vals_type), POINTER :: force_mixing_section, &
|
||||
link_section, qm_kind_section
|
||||
|
||||
!
|
||||
! create new qmmm sections for core and extended
|
||||
|
||||
NULLIFY(qmmm_core_section, qmmm_extended_section)
|
||||
CALL section_vals_duplicate(qmmm_section, qmmm_core_section, error=error)
|
||||
CALL section_vals_duplicate(qmmm_section, qmmm_extended_section, error=error)
|
||||
|
||||
! remove LINKs (specified by user for core) from extended
|
||||
link_section => section_vals_get_subs_vals(qmmm_extended_section, "LINK", can_return_null=.TRUE., error=error)
|
||||
IF (ASSOCIATED(link_section)) THEN
|
||||
CALL section_vals_remove_values(link_section,error)
|
||||
END IF
|
||||
|
||||
force_mixing_section => section_vals_get_subs_vals(qmmm_section,"FORCE_MIXING",error=error)
|
||||
|
||||
! get QM_KIND_ELEMENT_MAPPING
|
||||
CALL section_vals_val_get(force_mixing_section,"QM_KIND_ELEMENT_MAPPING",n_rep_val=n_elements,error=error)
|
||||
ALLOCATE(elem_mapping(2,n_elements))
|
||||
DO ielem=1,n_elements
|
||||
CALL section_vals_val_get(force_mixing_section,"QM_KIND_ELEMENT_MAPPING",i_rep_val=ielem,c_vals=elem_mapping_entry,&
|
||||
error=error)
|
||||
elem_mapping(1:2,ielem) = elem_mapping_entry(1:2)
|
||||
END DO
|
||||
|
||||
! get CUR_INDICES, CUR_LABELS
|
||||
CALL get_force_mixing_indices(force_mixing_section, cur_indices, cur_labels, error=error)
|
||||
CALL cp_assert(SIZE(cur_indices) > 0,cp_failure_level,cp_assertion_failed,routineP,&
|
||||
"cur_indices is empty, found no QM atoms"//&
|
||||
CPSourceFileRef)
|
||||
|
||||
! loop through elements and atoms, and set up new QM_KIND sections
|
||||
particles => subsys%particles%els
|
||||
|
||||
qm_kind_section => section_vals_get_subs_vals3(qmmm_section,"QM_KIND",error=error)
|
||||
CALL section_vals_get(qm_kind_section,n_repetition=i_rep_section_core,error=error)
|
||||
CALL cp_assert(i_rep_section_core > 0,cp_failure_level,cp_assertion_failed,routineP,&
|
||||
"Force-mixing QM didn't find any QM_KIND sections, so no core specified!"//&
|
||||
CPSourceFileRef)
|
||||
i_rep_section_extended = i_rep_section_core
|
||||
DO ielem=1,n_elements
|
||||
new_element_core = .TRUE.
|
||||
new_element_extended = .TRUE.
|
||||
DO ip=1, SIZE(cur_indices) ! particles with label
|
||||
|
||||
! extended
|
||||
! if current particle is some sort of QM atom, and not in core list (those the user
|
||||
! gave explicit QM_KIND sections for, and not a termination atom, need to make a
|
||||
! QM_KIND section for it
|
||||
IF (cur_labels(ip) > force_mixing_label_none .AND. cur_labels(ip) < force_mixing_label_QM_core_list .AND. &
|
||||
cur_labels(ip) /= force_mixing_label_termination) THEN
|
||||
IF (TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol) == TRIM(elem_mapping(1,ielem))) THEN
|
||||
qm_kind_section => section_vals_get_subs_vals3(qmmm_extended_section,"QM_KIND",error=error)
|
||||
IF (new_element_extended) THEN ! add new QM_KIND section for this element
|
||||
i_rep_section_extended = i_rep_section_extended + 1
|
||||
CALL section_vals_add_values(qm_kind_section,error=error)
|
||||
CALL section_vals_val_set(qm_kind_section,"_SECTION_PARAMETERS_",i_rep_section=i_rep_section_extended, &
|
||||
c_val=elem_mapping(2,ielem),error=error)
|
||||
i_rep_val_extended = 0
|
||||
new_element_extended = .FALSE.
|
||||
ENDIF
|
||||
i_rep_val_extended = i_rep_val_extended + 1
|
||||
CALL section_vals_val_set(qm_kind_section,"MM_INDEX",i_rep_section=i_rep_section_extended, &
|
||||
i_rep_val=i_rep_val_extended, i_val=cur_indices(ip),error=error)
|
||||
ENDIF ! element name matches current element
|
||||
ENDIF ! is a non-termination QM atom
|
||||
|
||||
! core
|
||||
! if current particle is a core QM atom, and not in core list (those the user
|
||||
! gave explicit QM_KIND sections for, need to make a QM_KIND section for it
|
||||
IF (cur_labels(ip) == force_mixing_label_QM_core) THEN
|
||||
IF (TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol) == TRIM(elem_mapping(1,ielem))) THEN
|
||||
qm_kind_section => section_vals_get_subs_vals3(qmmm_core_section,"QM_KIND",error=error)
|
||||
IF (new_element_core) THEN ! add new QM_KIND section for this element
|
||||
i_rep_section_core = i_rep_section_core + 1
|
||||
CALL section_vals_add_values(qm_kind_section,error=error)
|
||||
CALL section_vals_val_set(qm_kind_section,"_SECTION_PARAMETERS_",i_rep_section=i_rep_section_core, &
|
||||
c_val=elem_mapping(2,ielem),error=error)
|
||||
i_rep_val_core = 0
|
||||
new_element_core = .FALSE.
|
||||
ENDIF
|
||||
i_rep_val_core = i_rep_val_core + 1
|
||||
CALL section_vals_val_set(qm_kind_section,"MM_INDEX",i_rep_section=i_rep_section_core, &
|
||||
i_rep_val=i_rep_val_core, i_val=cur_indices(ip),error=error)
|
||||
ENDIF ! element name matches current element
|
||||
ENDIF ! is a non-termination QM atom
|
||||
|
||||
END DO ! atom index ip
|
||||
END DO ! element index ielem
|
||||
|
||||
![NB] check
|
||||
DEALLOCATE(elem_mapping, cur_indices, cur_labels)
|
||||
|
||||
END SUBROUTINE setup_force_mixing_qmmm_sections
|
||||
|
||||
SUBROUTINE get_force_mixing_indices(force_mixing_section, indices, labels, error)
|
||||
TYPE(section_vals_type), POINTER :: force_mixing_section
|
||||
INTEGER, POINTER :: indices(:), labels(:)
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'get_force_mixing_indices', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: i_rep_val, n_indices, &
|
||||
n_labels, n_reps
|
||||
INTEGER, POINTER :: indices_entry(:), &
|
||||
labels_entry(:)
|
||||
LOGICAL :: explicit
|
||||
TYPE(section_vals_type), POINTER :: restart_section
|
||||
|
||||
NULLIFY(indices, labels)
|
||||
restart_section => section_vals_get_subs_vals(force_mixing_section, "RESTART_INFO", error=error)
|
||||
CALL section_vals_get(restart_section,explicit=explicit,error=error)
|
||||
IF (.NOT. explicit) THEN ! no old indices, labels, return empty arrays
|
||||
ALLOCATE(indices(0))
|
||||
ALLOCATE(labels(0))
|
||||
RETURN
|
||||
ENDIF
|
||||
|
||||
![NB] maybe switch to reallocatable array
|
||||
CALL section_vals_val_get(restart_section,"INDICES",n_rep_val=n_reps,error=error)
|
||||
n_indices=0
|
||||
DO i_rep_val = 1,n_reps
|
||||
CALL section_vals_val_get(restart_section,"INDICES",&
|
||||
i_rep_val=i_rep_val,i_vals=indices_entry,error=error)
|
||||
n_indices = n_indices + SIZE(indices_entry)
|
||||
END DO
|
||||
ALLOCATE(indices(n_indices))
|
||||
n_indices=0
|
||||
DO i_rep_val = 1,n_reps
|
||||
CALL section_vals_val_get(restart_section,"INDICES",&
|
||||
i_rep_val=i_rep_val,i_vals=indices_entry,error=error)
|
||||
indices(n_indices+1:n_indices+SIZE(indices_entry)) = indices_entry
|
||||
n_indices = n_indices + SIZE(indices_entry)
|
||||
END DO
|
||||
|
||||
CALL section_vals_val_get(restart_section,"LABELS",n_rep_val=n_reps,error=error)
|
||||
n_labels=0
|
||||
DO i_rep_val = 1,n_reps
|
||||
CALL section_vals_val_get(restart_section,"LABELS",&
|
||||
i_rep_val=i_rep_val,i_vals=labels_entry,error=error)
|
||||
n_labels = n_labels + SIZE(labels_entry)
|
||||
END DO
|
||||
ALLOCATE(labels(n_labels))
|
||||
n_labels=0
|
||||
DO i_rep_val = 1,n_reps
|
||||
CALL section_vals_val_get(restart_section,"LABELS",&
|
||||
i_rep_val=i_rep_val,i_vals=labels_entry,error=error)
|
||||
labels(n_labels+1:n_labels+SIZE(labels_entry)) = labels_entry
|
||||
n_labels = n_labels + SIZE(labels_entry)
|
||||
END DO
|
||||
|
||||
CALL cp_assert(n_indices == n_labels, cp_failure_level,cp_assertion_failed,routineP,&
|
||||
"got unequal numbers of force_mixing indices and labels!"//&
|
||||
CPSourceFileRef)
|
||||
END SUBROUTINE get_force_mixing_indices
|
||||
|
||||
END MODULE qmmm_force_mixing
|
||||
210
src/qmmm_main.F
210
src/qmmm_main.F
|
|
@ -23,6 +23,7 @@ MODULE qmmm_main
|
|||
USE cp_para_types, ONLY: cp_para_env_type
|
||||
USE cp_subsys_methods, ONLY: create_small_subsys
|
||||
USE cp_subsys_types, ONLY: cp_subsys_release,&
|
||||
cp_subsys_retain,&
|
||||
cp_subsys_type
|
||||
USE f77_blas
|
||||
USE fist_environment_types, ONLY: fist_env_get,&
|
||||
|
|
@ -32,6 +33,7 @@ MODULE qmmm_main
|
|||
USE force_env_methods, ONLY: force_env_create
|
||||
USE force_env_types, ONLY: force_env_get,&
|
||||
force_env_p_type,&
|
||||
force_env_release,&
|
||||
force_env_type
|
||||
USE global_types, ONLY: global_environment_type
|
||||
USE header, ONLY: qmmm_header
|
||||
|
|
@ -39,13 +41,18 @@ MODULE qmmm_main
|
|||
do_qmmm,&
|
||||
do_qmmm_none,&
|
||||
do_qs
|
||||
USE input_cp2k_restarts, ONLY: update_input
|
||||
USE input_section_types, ONLY: section_vals_get,&
|
||||
section_vals_get_subs_vals,&
|
||||
section_vals_release,&
|
||||
section_vals_type,&
|
||||
section_vals_val_get,&
|
||||
section_vals_val_set
|
||||
USE kinds, ONLY: default_string_length,&
|
||||
dp
|
||||
USE md_environment_types, ONLY: md_environment_type
|
||||
USE qmmm_force_mixing, ONLY: setup_force_mixing_qmmm_sections,&
|
||||
update_force_mixing_labels
|
||||
USE qmmm_init, ONLY: &
|
||||
assign_mm_charges_and_radius, move_or_add_atoms, &
|
||||
print_image_charge_info, print_qmmm_charges, print_qmmm_links, &
|
||||
|
|
@ -55,15 +62,19 @@ MODULE qmmm_main
|
|||
USE qmmm_links_methods, ONLY: qmmm_link_Imomm_coord
|
||||
USE qmmm_pw_grid, ONLY: qmmm_pw_grid_init
|
||||
USE qmmm_types, ONLY: &
|
||||
add_set_release, add_set_type, qmmm_env_mm_create, &
|
||||
add_set_release, add_set_type, fist_subsys, force_mixing_core_subsys, &
|
||||
force_mixing_extended_subsys, primary_subsys, qmmm_env_mm_create, &
|
||||
qmmm_env_mm_release, qmmm_env_mm_type, qmmm_env_qm_create, &
|
||||
qmmm_env_qm_release, qmmm_env_qm_type, qmmm_links_type
|
||||
qmmm_env_qm_release, qmmm_env_qm_retain, qmmm_env_qm_type, &
|
||||
qmmm_links_type, qs_subsys
|
||||
USE qs_environment_types, ONLY: qs_environment_type,&
|
||||
set_qs_env
|
||||
USE qs_main, ONLY: quickstep_create_force_env
|
||||
USE termination, ONLY: stop_program
|
||||
USE timings, ONLY: timeset,&
|
||||
timestop
|
||||
USE virial_types, ONLY: virial_release,&
|
||||
virial_retain
|
||||
#include "cp_common_uses.h"
|
||||
|
||||
IMPLICIT NONE
|
||||
|
|
@ -72,10 +83,182 @@ MODULE qmmm_main
|
|||
LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE.
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_main'
|
||||
|
||||
PUBLIC :: qmmm_create_force_env
|
||||
PUBLIC :: qmmm_create_force_env, qmmm_update_force_mixing_env
|
||||
|
||||
|
||||
!***
|
||||
CONTAINS
|
||||
|
||||
! *****************************************************************************
|
||||
!> \param error variable to control error logging, stopping,...
|
||||
!> see module cp_error_handling
|
||||
!> \par History
|
||||
!> 02.2012 created [noam]
|
||||
!> \author Noam Bernstein
|
||||
! *****************************************************************************
|
||||
SUBROUTINE qmmm_create_force_env(force_env, root_section, para_env, globenv,&
|
||||
force_env_section, subsys_section, use_motion_section, dummy_force_env, error)
|
||||
TYPE(force_env_type), POINTER :: force_env
|
||||
TYPE(section_vals_type), POINTER :: root_section
|
||||
TYPE(cp_para_env_type), POINTER :: para_env
|
||||
TYPE(global_environment_type), POINTER :: globenv
|
||||
TYPE(section_vals_type), POINTER :: force_env_section, &
|
||||
subsys_section
|
||||
LOGICAL, INTENT(IN) :: use_motion_section
|
||||
TYPE(force_env_type), OPTIONAL, POINTER :: dummy_force_env
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_create_force_env', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: stat
|
||||
LOGICAL :: failure, force_mixing_active
|
||||
TYPE(cp_subsys_type), POINTER :: subsys
|
||||
TYPE(force_env_p_type), DIMENSION(:), &
|
||||
POINTER :: sub_force_env
|
||||
TYPE(force_env_type), POINTER :: use_dummy_force_env
|
||||
TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm
|
||||
TYPE(section_vals_type), POINTER :: qmmm_core_section, &
|
||||
qmmm_extended_section, &
|
||||
qmmm_force_mixing, &
|
||||
qmmm_section
|
||||
|
||||
qmmm_section => section_vals_get_subs_vals(force_env_section,"QMMM",error=error)
|
||||
qmmm_force_mixing => section_vals_get_subs_vals(qmmm_section,"FORCE_MIXING",error=error)
|
||||
CALL section_vals_get(qmmm_force_mixing,explicit=force_mixing_active,error=error)
|
||||
|
||||
IF (force_mixing_active) THEN
|
||||
IF (PRESENT(dummy_force_env)) THEN
|
||||
use_dummy_force_env => force_env
|
||||
ELSE
|
||||
CALL qmmm_create_force_env_low(use_dummy_force_env, &
|
||||
root_section, para_env, globenv,&
|
||||
force_env_section, qmmm_section, subsys_section, use_motion_section, error)
|
||||
ENDIF
|
||||
CALL force_env_get(use_dummy_force_env, subsys=subsys, error=error)
|
||||
|
||||
![NB]
|
||||
! check that nothing accesses via root or force_env sections
|
||||
ALLOCATE(sub_force_env(2),stat=stat)
|
||||
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
|
||||
|
||||
CALL update_force_mixing_labels(use_dummy_force_env, subsys, qmmm_section, error=error)
|
||||
|
||||
! using CUR_INDICES and CUR_LABELS, create appropriate QM_KIND sections for two QM/MM calculations
|
||||
CALL setup_force_mixing_qmmm_sections(subsys, qmmm_section, qmmm_core_section, qmmm_extended_section,error=error)
|
||||
|
||||
NULLIFY(sub_force_env(force_mixing_core_subsys)%force_env)
|
||||
CALL qmmm_create_force_env_low(sub_force_env(force_mixing_core_subsys)%force_env, &
|
||||
root_section, para_env, globenv,&
|
||||
force_env_section, qmmm_core_section, subsys_section, use_motion_section, error)
|
||||
NULLIFY(sub_force_env(force_mixing_extended_subsys)%force_env)
|
||||
CALL qmmm_create_force_env_low(sub_force_env(force_mixing_extended_subsys)%force_env, &
|
||||
root_section, para_env, globenv,&
|
||||
force_env_section, qmmm_extended_section, subsys_section, use_motion_section, error)
|
||||
|
||||
CALL section_vals_release(qmmm_core_section, error=error)
|
||||
CALL section_vals_release(qmmm_extended_section, error=error)
|
||||
|
||||
IF (.NOT. PRESENT(dummy_force_env)) THEN ! allocated it above, now release
|
||||
CALL force_env_release(use_dummy_force_env, error=error)
|
||||
ENDIF
|
||||
ELSE
|
||||
ALLOCATE(sub_force_env(1),stat=stat)
|
||||
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
|
||||
CALL qmmm_create_force_env_low(sub_force_env(1)%force_env, &
|
||||
root_section, para_env, globenv,&
|
||||
force_env_section, qmmm_section, subsys_section, use_motion_section, error)
|
||||
ENDIF
|
||||
|
||||
CALL force_env_get(sub_force_env(primary_subsys)%force_env,qmmm_env=qmmm_env_qm,error=error)
|
||||
CALL force_env_create(force_env, root_section, sub_force_env=sub_force_env,&
|
||||
qmmm_env=qmmm_env_qm, &
|
||||
para_env=para_env, globenv=globenv, force_env_section= force_env_section,&
|
||||
error=error)
|
||||
|
||||
END SUBROUTINE qmmm_create_force_env
|
||||
|
||||
|
||||
SUBROUTINE qmmm_update_force_mixing_env(md_env, force_env, root_section, error)
|
||||
TYPE(md_environment_type), POINTER :: md_env
|
||||
TYPE(force_env_type), POINTER :: force_env
|
||||
TYPE(section_vals_type), POINTER :: root_section
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_update_force_mixing_env', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: stat
|
||||
LOGICAL :: failure, force_mixing_active, &
|
||||
labels_changed
|
||||
TYPE(cp_subsys_type), POINTER :: subsys
|
||||
TYPE(force_env_p_type), DIMENSION(:), &
|
||||
POINTER :: sub_force_env
|
||||
TYPE(section_vals_type), POINTER :: qmmm_core_section, &
|
||||
qmmm_extended_Section, qmmm_force_mixing, qmmm_section, subsys_section
|
||||
|
||||
qmmm_section => section_vals_get_subs_vals(force_env%force_env_section,"QMMM",error=error)
|
||||
qmmm_force_mixing => section_vals_get_subs_vals(qmmm_section,"FORCE_MIXING",error=error)
|
||||
CALL section_vals_get(qmmm_force_mixing,explicit=force_mixing_active,error=error)
|
||||
|
||||
IF (force_mixing_active) THEN
|
||||
|
||||
CALL force_env_get(force_env, subsys=subsys, error=error)
|
||||
CALL update_force_mixing_labels(force_env, subsys, qmmm_section, labels_changed=labels_changed, error=error)
|
||||
IF (.NOT. labels_changed) THEN
|
||||
RETURN
|
||||
ENDIF
|
||||
CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,&
|
||||
"Adaptive force-mixing labels changed, rebuilding QM/MM calculations! "//&
|
||||
CPSourceFileRef,&
|
||||
only_ionode=.TRUE.)
|
||||
|
||||
CALL update_input(md_env,force_env,root_section,error=error)
|
||||
|
||||
ALLOCATE(sub_force_env(2),stat=stat)
|
||||
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
|
||||
|
||||
! using CUR_INDICES and CUR_LABELS, create appropriate QM_KIND sections for two QM/MM calculations
|
||||
CALL setup_force_mixing_qmmm_sections(subsys, qmmm_section, qmmm_core_section, qmmm_extended_section,error=error)
|
||||
|
||||
subsys_section => section_vals_get_subs_vals(force_env%force_env_section,"SUBSYS",error=error)
|
||||
NULLIFY(sub_force_env(force_mixing_core_subsys)%force_env)
|
||||
![ADAPT] no sure about use_motion_section
|
||||
CALL qmmm_create_force_env_low(sub_force_env(force_mixing_core_subsys)%force_env, &
|
||||
force_env%root_section, force_env%para_env, force_env%globenv,&
|
||||
force_env%force_env_section, qmmm_core_section, subsys_section, use_motion_section=.TRUE., error=error)
|
||||
NULLIFY(sub_force_env(force_mixing_extended_subsys)%force_env)
|
||||
CALL qmmm_create_force_env_low(sub_force_env(force_mixing_extended_subsys)%force_env, &
|
||||
force_env%root_section, force_env%para_env, force_env%globenv,&
|
||||
force_env%force_env_section, qmmm_extended_section, subsys_section, use_motion_section=.TRUE., error=error)
|
||||
|
||||
CALL section_vals_release(qmmm_core_section, error=error)
|
||||
CALL section_vals_release(qmmm_extended_section, error=error)
|
||||
|
||||
! do something with sub_force_envs
|
||||
|
||||
! duplicate releases in force_env_release
|
||||
CALL force_env_release(force_env%sub_force_env(force_mixing_core_subsys)%force_env,error=error)
|
||||
CALL force_env_release(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,error=error)
|
||||
CALL qmmm_env_qm_release(force_env%qmmm_env,error=error)
|
||||
CALL virial_release(force_env%virial,error=error)
|
||||
CALL cp_subsys_release(force_env%subsys,error=error)
|
||||
|
||||
DEALLOCATE(force_env%sub_force_env)
|
||||
force_env%sub_force_env => sub_force_env
|
||||
|
||||
! duplicate things that are pointed to in force_env_create when doing qmmm
|
||||
force_env%virial => sub_force_env(primary_subsys)%force_env%virial
|
||||
CALL virial_retain(force_env%virial, error=error)
|
||||
force_env%qmmm_env => sub_force_env(primary_subsys)%force_env%qmmm_env
|
||||
CALL qmmm_env_qm_retain(force_env%qmmm_env, error=error)
|
||||
force_env%subsys => sub_force_env(primary_subsys)%force_env%subsys
|
||||
CALL cp_subsys_retain(force_env%subsys, error=error)
|
||||
|
||||
ENDIF
|
||||
|
||||
END SUBROUTINE qmmm_update_force_mixing_env
|
||||
|
||||
! *****************************************************************************
|
||||
!> \param error variable to control error logging, stopping,...
|
||||
!> see module cp_error_handling
|
||||
|
|
@ -83,18 +266,18 @@ CONTAINS
|
|||
!> 05.2004 created [fawzi]
|
||||
!> \author Fawzi Mohamed
|
||||
! *****************************************************************************
|
||||
SUBROUTINE qmmm_create_force_env(force_env, root_section, para_env, globenv,&
|
||||
force_env_section, subsys_section, use_motion_section, error)
|
||||
SUBROUTINE qmmm_create_force_env_low(force_env, root_section, para_env, globenv,&
|
||||
force_env_section, qmmm_section, subsys_section, use_motion_section, error)
|
||||
TYPE(force_env_type), POINTER :: force_env
|
||||
TYPE(section_vals_type), POINTER :: root_section
|
||||
TYPE(cp_para_env_type), POINTER :: para_env
|
||||
TYPE(global_environment_type), POINTER :: globenv
|
||||
TYPE(section_vals_type), POINTER :: force_env_section, &
|
||||
subsys_section
|
||||
qmmm_section, subsys_section
|
||||
LOGICAL, INTENT(IN) :: use_motion_section
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_create_force_env', &
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_create_force_env_low', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
CHARACTER(len=default_string_length), &
|
||||
|
|
@ -126,7 +309,7 @@ CONTAINS
|
|||
TYPE(qmmm_links_type), POINTER :: qmmm_links
|
||||
TYPE(qs_environment_type), POINTER :: qs_env
|
||||
TYPE(section_vals_type), POINTER :: print_gen, print_section, &
|
||||
qmmm_periodic, qmmm_section
|
||||
qmmm_periodic
|
||||
|
||||
CALL timeset(routineN,handle)
|
||||
failure=.FALSE.
|
||||
|
|
@ -146,7 +329,6 @@ CONTAINS
|
|||
subsys_section => section_vals_get_subs_vals(force_env_section,"SUBSYS",&
|
||||
error=error)
|
||||
END IF
|
||||
qmmm_section => section_vals_get_subs_vals(force_env_section,"QMMM",error=error)
|
||||
qmmm_periodic => section_vals_get_subs_vals(qmmm_section,"PERIODIC",error=error)
|
||||
print_section => section_vals_get_subs_vals(qmmm_section,"PRINT",error=error)
|
||||
print_gen => section_vals_get_subs_vals(print_section,"PROGRAM_RUN_INFO",error=error)
|
||||
|
|
@ -209,12 +391,12 @@ CPSourceFileRef,&
|
|||
|
||||
! First Initialize Fist...
|
||||
CALL section_vals_val_set(force_env_section,"METHOD",i_val=do_fist,error=error)
|
||||
CALL fist_create_force_env(sub_force_env(1)%force_env, &
|
||||
CALL fist_create_force_env(sub_force_env(fist_subsys)%force_env, &
|
||||
root_section, para_env, globenv, qmmm=.TRUE., qmmm_env=qmmm_env_mm,&
|
||||
force_env_section=force_env_section, subsys_section=subsys_section,&
|
||||
use_motion_section=use_motion_section,error=error)
|
||||
|
||||
CALL force_env_get(sub_force_env(1)%force_env,subsys=subsys_mm,&
|
||||
CALL force_env_get(sub_force_env(fist_subsys)%force_env,subsys=subsys_mm,&
|
||||
cell=mm_cell,fist_env=fist_env,error=error)
|
||||
|
||||
! Set up QM/MM Options
|
||||
|
|
@ -328,13 +510,13 @@ CPSourceFileRef,&
|
|||
force_env_section=force_env_section, subsys_section=subsys_section,error=error)
|
||||
IF (qmmm_link_imomm) CALL qmmm_link_Imomm_coord(qmmm_links, subsys_qm%particles%els,&
|
||||
qm_atom_index, error)
|
||||
CALL quickstep_create_force_env(sub_force_env(2)%force_env,root_section,para_env,&
|
||||
CALL quickstep_create_force_env(sub_force_env(qs_subsys)%force_env,root_section,para_env,&
|
||||
globenv, subsys=subsys_qm, cell=qm_cell_small, qmmm=.TRUE.,qmmm_env_qm=qmmm_env_qm,&
|
||||
qmmm_periodic=qmmm_env_qm%periodic, force_env_section=force_env_section,&
|
||||
subsys_section=subsys_section,use_motion_section=use_motion_section,&
|
||||
error=error)
|
||||
CALL cp_subsys_release(subsys_qm,error=error)
|
||||
CALL force_env_get(sub_force_env(2)%force_env,qs_env=qs_env,error=error)
|
||||
CALL force_env_get(sub_force_env(qs_subsys)%force_env,qs_env=qs_env,error=error)
|
||||
IF (qmmm_env_qm%periodic) THEN
|
||||
IF (.NOT.ASSOCIATED(super_cell)) THEN
|
||||
ALLOCATE (super_cell,STAT=stat)
|
||||
|
|
@ -414,6 +596,6 @@ CPSourceFileRef,&
|
|||
END IF
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE qmmm_create_force_env
|
||||
END SUBROUTINE qmmm_create_force_env_low
|
||||
|
||||
END MODULE qmmm_main
|
||||
|
|
|
|||
|
|
@ -35,8 +35,17 @@ MODULE qmmm_types
|
|||
LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE.
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_types'
|
||||
INTEGER, SAVE, PRIVATE :: last_qmmm_env_id_nr=0
|
||||
INTEGER, PARAMETER, PUBLIC :: primary_subsys = 1
|
||||
INTEGER, PARAMETER, PUBLIC :: fist_subsys = 1,&
|
||||
qs_subsys = 2
|
||||
INTEGER, PARAMETER, PUBLIC :: force_mixing_extended_subsys = 1,&
|
||||
force_mixing_core_subsys = 2
|
||||
INTEGER, PARAMETER, PUBLIC :: force_mixing_label_none = -1,&
|
||||
force_mixing_label_QM_core_list = 10,&
|
||||
force_mixing_label_QM_core = 9,&
|
||||
force_mixing_label_QM_extended = 8,&
|
||||
force_mixing_label_buffer = 7,&
|
||||
force_mixing_label_termination = 6
|
||||
|
||||
PUBLIC :: qmmm_env_qm_type, gridlevel_info_type, qmmm_pot_type, qmmm_pot_p_type
|
||||
PUBLIC :: qmmm_env_qm_retain, qmmm_env_qm_release, qmmm_env_qm_create
|
||||
|
|
|
|||
139
src/qmmm_util.F
139
src/qmmm_util.F
|
|
@ -9,13 +9,15 @@
|
|||
!> \author Teodoro Laino
|
||||
! *****************************************************************************
|
||||
MODULE qmmm_util
|
||||
USE cell_types, ONLY: cell_type
|
||||
USE cell_types, ONLY: cell_copy,&
|
||||
cell_type
|
||||
USE cp_subsys_types, ONLY: cp_subsys_type
|
||||
USE f77_blas
|
||||
USE force_env_types, ONLY: force_env_get,&
|
||||
force_env_type,&
|
||||
use_qmmm
|
||||
USE input_constants, ONLY: do_qmmm_wall_quadratic,&
|
||||
USE input_constants, ONLY: do_qmmm_wall_none,&
|
||||
do_qmmm_wall_quadratic,&
|
||||
do_qmmm_wall_reflective
|
||||
USE input_section_types, ONLY: section_vals_get,&
|
||||
section_vals_get_subs_vals,&
|
||||
|
|
@ -23,8 +25,12 @@ MODULE qmmm_util
|
|||
section_vals_val_get
|
||||
USE kinds, ONLY: dp
|
||||
USE particle_types, ONLY: particle_type,&
|
||||
write_fist_particle_coordinates,&
|
||||
write_qs_particle_coordinates
|
||||
USE qmmm_types, ONLY: fist_subsys,&
|
||||
force_mixing_core_subsys,&
|
||||
force_mixing_extended_subsys,&
|
||||
primary_subsys,&
|
||||
qs_subsys
|
||||
USE qs_energy_types, ONLY: qs_energy_type
|
||||
USE qs_environment_types, ONLY: get_qs_env
|
||||
|
|
@ -33,11 +39,12 @@ MODULE qmmm_util
|
|||
IMPLICIT NONE
|
||||
PRIVATE
|
||||
|
||||
LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE.
|
||||
LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.FALSE.
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_util'
|
||||
PUBLIC :: apply_qmmm_walls_reflective,&
|
||||
apply_qmmm_walls,&
|
||||
apply_qmmm_translate,&
|
||||
apply_qmmm_translate_low,&
|
||||
spherical_cutoff_factor
|
||||
|
||||
CONTAINS
|
||||
|
|
@ -119,6 +126,8 @@ CONTAINS
|
|||
CALL section_vals_val_get(walls_section,"TYPE",i_val=iwall_type,error=error)
|
||||
skin(:) = list(:)
|
||||
ELSE
|
||||
![NB]
|
||||
iwall_type=do_qmmm_wall_reflective
|
||||
skin(:) = 0.0_dp
|
||||
END IF
|
||||
CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure)
|
||||
|
|
@ -126,9 +135,16 @@ CONTAINS
|
|||
CPPrecondition(ASSOCIATED(force_env%qmmm_env),cp_failure_level,routineP,error,failure)
|
||||
CPPrecondition(force_env%qmmm_env%ref_count>0,cp_failure_level,routineP,error,failure)
|
||||
|
||||
CALL force_env_get(force_env%sub_force_env(fist_subsys)%force_env,&
|
||||
CALL cp_assert(iwall_type == do_qmmm_wall_none .OR. SIZE(force_env%sub_force_env) == 1, &
|
||||
cp_warning_level,cp_assertion_failed,routineP,&
|
||||
"Reflective walls for QM/MM are not implemented (or useful), when "//&
|
||||
"top level force_env has two or more sub_force_env allocated, which only happens when "//&
|
||||
"force mixing is active. Skipping!"//CPSourceFileRef)
|
||||
IF (SIZE(force_env%sub_force_env) /= 1) RETURN
|
||||
|
||||
CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(fist_subsys)%force_env,&
|
||||
cell=mm_cell,subsys=subsys_mm,error=error)
|
||||
CALL force_env_get(force_env%sub_force_env(qs_subsys)%force_env,&
|
||||
CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%force_env,&
|
||||
cell=qm_cell,subsys=subsys_qm,error=error)
|
||||
qm_atom_index => force_env%qmmm_env%qm_atom_index
|
||||
CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure)
|
||||
|
|
@ -228,10 +244,19 @@ CONTAINS
|
|||
CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure)
|
||||
CPPrecondition(ASSOCIATED(force_env%qmmm_env),cp_failure_level,routineP,error,failure)
|
||||
CPPrecondition(force_env%qmmm_env%ref_count>0,cp_failure_level,routineP,error,failure)
|
||||
CALL force_env_get(force_env%sub_force_env(fist_subsys)%force_env,&
|
||||
|
||||
CALL cp_assert(SIZE(force_env%sub_force_env) == 1, &
|
||||
cp_warning_level,cp_assertion_failed,routineP,&
|
||||
"Quadratic walls for QM/MM are not implemented (or useful), when "//&
|
||||
"top level force_env has two or more sub_force_env allocated, which only happens when "//&
|
||||
"force mixing is active. Skipping!"//CPSourceFileRef)
|
||||
IF (SIZE(force_env%sub_force_env) /= 1) RETURN
|
||||
|
||||
CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(fist_subsys)%force_env,&
|
||||
cell=mm_cell,subsys=subsys_mm,error=error)
|
||||
CALL force_env_get(force_env%sub_force_env(qs_subsys)%force_env,&
|
||||
CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%force_env,&
|
||||
cell=qm_cell,subsys=subsys_qm,error=error)
|
||||
|
||||
qm_atom_index => force_env%qmmm_env%qm_atom_index
|
||||
CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure)
|
||||
|
||||
|
|
@ -294,11 +319,11 @@ CONTAINS
|
|||
ENDIF
|
||||
ENDDO
|
||||
|
||||
CALL get_qs_env(qs_env=force_env%sub_force_env(qs_subsys)%&
|
||||
CALL get_qs_env(qs_env=force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%&
|
||||
force_env%qs_env, energy=energy,error=error)
|
||||
energy%total = energy%total + wallenergy
|
||||
|
||||
END SUBROUTINE
|
||||
END SUBROUTINE apply_qmmm_walls_quadratic
|
||||
|
||||
! *****************************************************************************
|
||||
!> \brief Apply translation to the full system in order to center the QM
|
||||
|
|
@ -316,9 +341,95 @@ CONTAINS
|
|||
CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_translate', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: ip
|
||||
INTEGER, DIMENSION(:), POINTER :: qm_atom_index
|
||||
LOGICAL :: failure
|
||||
REAL(dp), POINTER :: charges(:)
|
||||
TYPE(cell_type), POINTER :: cell_core, cell_extended
|
||||
TYPE(cp_subsys_type), POINTER :: subsys_core, subsys_extended, &
|
||||
subsys_fist, subsys_qs
|
||||
TYPE(particle_type), DIMENSION(:), &
|
||||
POINTER :: particles_core, &
|
||||
particles_extended, &
|
||||
particles_mm, particles_qm
|
||||
TYPE(section_vals_type), POINTER :: subsys_section
|
||||
|
||||
SELECT CASE(force_env%in_use)
|
||||
|
||||
CASE(use_qmmm)
|
||||
IF (SIZE(force_env%sub_force_env) == 1) THEN
|
||||
CALL apply_qmmm_translate_low(force_env%sub_force_env(1)%force_env,error)
|
||||
ELSE IF (SIZE(force_env%sub_force_env) == 2) THEN
|
||||
! want to center extended, and make core consistent with that
|
||||
CALL apply_qmmm_translate_low(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,error)
|
||||
|
||||
! translate core fist particles
|
||||
CALL force_env_get(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,&
|
||||
subsys=subsys_extended,cell=cell_extended,error=error)
|
||||
CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env,&
|
||||
subsys=subsys_core,cell=cell_core,error=error)
|
||||
particles_extended => subsys_extended%particles%els
|
||||
particles_core => subsys_core%particles%els
|
||||
DO ip=1,SIZE(particles_extended)
|
||||
particles_core(ip)%r = particles_extended(ip)%r
|
||||
END DO
|
||||
CALL cell_copy(cell_extended, cell_core, error)
|
||||
|
||||
! also update core QM particles (like in apply_qmmm_translate_low)
|
||||
CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(fist_subsys)%force_env,&
|
||||
subsys=subsys_fist,error=error)
|
||||
CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(qs_subsys)%force_env,&
|
||||
subsys=subsys_qs,error=error)
|
||||
qm_atom_index => force_env%sub_force_env(force_mixing_core_subsys)%force_env%qmmm_env%qm_atom_index
|
||||
particles_qm => subsys_qs%particles%els
|
||||
particles_mm => subsys_fist%particles%els
|
||||
CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure)
|
||||
DO ip=1,SIZE(qm_atom_index)
|
||||
particles_qm(ip)%r=particles_mm(qm_atom_index(ip))%r
|
||||
END DO
|
||||
|
||||
IF (debug_this_module) THEN
|
||||
CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(qs_subsys)%force_env,&
|
||||
subsys=subsys_core,cell=cell_core,error=error)
|
||||
particles_core => subsys_core%particles%els
|
||||
subsys_section => section_vals_get_subs_vals(force_env%force_env_section, &
|
||||
"SUBSYS",error=error)
|
||||
CALL write_qs_particle_coordinates(particles_core,subsys_section,"QM/MM core calc first QM, then MM (0 charges)",error)
|
||||
CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(fist_subsys)%force_env,&
|
||||
subsys=subsys_core,cell=cell_core,error=error)
|
||||
particles_core => subsys_core%particles%els
|
||||
ALLOCATE(charges(SIZE(particles_core)))
|
||||
charges = 0.0_dp
|
||||
CALL write_fist_particle_coordinates(particles_core,subsys_section,charges,error)
|
||||
DEALLOCATE(charges)
|
||||
ENDIF
|
||||
|
||||
ELSE
|
||||
CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure)
|
||||
ENDIF
|
||||
END SELECT
|
||||
END SUBROUTINE apply_qmmm_translate
|
||||
|
||||
! *****************************************************************************
|
||||
!> \brief Apply translation to the full system in order to center the QM
|
||||
!> system into the QM box
|
||||
!> \param error variable to control error logging, stopping,...
|
||||
!> see module cp_error_handling
|
||||
!> \par History
|
||||
!> 08.2007 created [tlaino] - Zurich University
|
||||
!> \author Teodoro Laino
|
||||
! *****************************************************************************
|
||||
SUBROUTINE apply_qmmm_translate_low(force_env,error)
|
||||
TYPE(force_env_type), POINTER :: force_env
|
||||
TYPE(cp_error_type), INTENT(inout) :: error
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_translate_low', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: ip, unit_nr
|
||||
INTEGER, DIMENSION(:), POINTER :: qm_atom_index
|
||||
LOGICAL :: failure
|
||||
REAL(DP), POINTER :: charges(:)
|
||||
REAL(KIND=dp), DIMENSION(3) :: max_coord, min_coord, transl_v
|
||||
TYPE(cell_type), POINTER :: mm_cell, qm_cell
|
||||
TYPE(cp_logger_type), POINTER :: logger
|
||||
|
|
@ -383,12 +494,18 @@ CONTAINS
|
|||
DO ip=1,SIZE(qm_atom_index)
|
||||
particles_qm(ip)%r=particles_mm(qm_atom_index(ip))%r
|
||||
END DO
|
||||
|
||||
subsys_section => section_vals_get_subs_vals(force_env%force_env_section, &
|
||||
"SUBSYS",error=error)
|
||||
CALL write_qs_particle_coordinates(particles_qm,subsys_section,"QM/MM",error)
|
||||
CALL write_qs_particle_coordinates(particles_qm,subsys_section,"QM/MM first QM, then MM (0 charges)",error)
|
||||
ALLOCATE(charges(SIZE(particles_mm)))
|
||||
charges = 0.0_dp
|
||||
CALL write_fist_particle_coordinates(particles_mm,subsys_section,charges,error)
|
||||
DEALLOCATE(charges)
|
||||
|
||||
END SELECT
|
||||
|
||||
END SUBROUTINE apply_qmmm_translate
|
||||
END SUBROUTINE apply_qmmm_translate_low
|
||||
|
||||
! *****************************************************************************
|
||||
!> \brief Computes a spherical cutoff factor for the QMMM interactions
|
||||
|
|
|
|||
122
tests/QMMM/QS/regtest-1/H2O-qmmm-none-force-mixing-1.inp
Normal file
122
tests/QMMM/QS/regtest-1/H2O-qmmm-none-force-mixing-1.inp
Normal file
|
|
@ -0,0 +1,122 @@
|
|||
&FORCE_EVAL
|
||||
METHOD QMMM
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME ../BASIS_SET
|
||||
POTENTIAL_FILE_NAME ../POTENTIAL
|
||||
&MGRID
|
||||
COMMENSURATE
|
||||
CUTOFF 50
|
||||
&END MGRID
|
||||
&QS
|
||||
&END QS
|
||||
&SCF
|
||||
SCF_GUESS atomic
|
||||
&END SCF
|
||||
&XC
|
||||
&XC_FUNCTIONAL pade
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&MM
|
||||
&FORCEFIELD
|
||||
parm_file_name ../sample_pot/5H2O.pot
|
||||
parmtype CHM
|
||||
&spline
|
||||
rcut_nb 5.2917720830
|
||||
&end
|
||||
&END FORCEFIELD
|
||||
&POISSON
|
||||
&EWALD
|
||||
EWALD_TYPE spme
|
||||
ALPHA .44
|
||||
GMAX 64
|
||||
O_SPLINE 6
|
||||
&END EWALD
|
||||
&END POISSON
|
||||
&END MM
|
||||
&QMMM
|
||||
# NOCENTER0
|
||||
NOCENTER
|
||||
MM_POTENTIAL_FILE_NAME ../MM_POTENTIAL
|
||||
&CELL
|
||||
ABC 17.0 7.0 7.0
|
||||
&END CELL
|
||||
ECOUPL GAUSS
|
||||
USE_GEEP_LIB 6
|
||||
NOCOMPATIBILITY
|
||||
|
||||
# core list
|
||||
&QM_KIND H
|
||||
MM_INDEX 2 3
|
||||
&END QM_KIND
|
||||
&QM_KIND O
|
||||
MM_INDEX 1
|
||||
&END QM_KIND
|
||||
|
||||
# Noam's new section
|
||||
&FORCE_MIXING
|
||||
R_CORE 2.5 3.5
|
||||
R_QM 2.5 3.5
|
||||
R_BUF 2.5 3.5
|
||||
QM_KIND_ELEMENT_MAPPING H H
|
||||
QM_KIND_ELEMENT_MAPPING O O
|
||||
# &RESTART_INFO
|
||||
# INDICES 1 2 3 4 5 6
|
||||
# LABELS 10 10 10 8 8 8
|
||||
# &END RESTART_INFO
|
||||
&END FORCE_MIXING
|
||||
#
|
||||
# QM_KINDS
|
||||
#
|
||||
&END QMMM
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 20.00 20.00 20.00
|
||||
&END CELL
|
||||
&KIND H
|
||||
BASIS_SET DZVP-GTH-BLYP
|
||||
POTENTIAL GTH-BLYP-q1
|
||||
&END KIND
|
||||
&KIND O
|
||||
BASIS_SET DZVP-GTH-BLYP
|
||||
POTENTIAL GTH-BLYP-q6
|
||||
&END KIND
|
||||
&TOPOLOGY
|
||||
CONN_FILE_NAME ../sample_psf/5H2O.psf
|
||||
CONNECTIVITY UPSF
|
||||
COORD_FILE_NAME ../sample_pdb/5H2O.pdb
|
||||
COORDINATE pdb
|
||||
&END TOPOLOGY
|
||||
&PRINT
|
||||
# &ATOMIC_COORDINATES
|
||||
# &END ATOMIC_COORDINATES
|
||||
# &CELL
|
||||
# &END CELL
|
||||
&END PRINT
|
||||
&END SUBSYS
|
||||
&PRINT
|
||||
&FORCES
|
||||
&END FORCES
|
||||
&END PRINT
|
||||
&END FORCE_EVAL
|
||||
&GLOBAL
|
||||
PROJECT H2O-qmmm-none-force_mixing-1
|
||||
PRINT_LEVEL MEDIUM
|
||||
RUN_TYPE MD
|
||||
&END GLOBAL
|
||||
&MOTION
|
||||
&MD
|
||||
STEPS 8
|
||||
ENSEMBLE REFTRAJ
|
||||
&REFTRAJ
|
||||
TRAJ_FILE_NAME ../sample_xyz/5H2O_adaptive.xyz
|
||||
EVAL_ENERGY_FORCES
|
||||
&END REFTRAJ
|
||||
&END MD
|
||||
|
||||
&PRINT
|
||||
&RESTART
|
||||
&END RESTART
|
||||
&END PRINT
|
||||
|
||||
&END MOTION
|
||||
|
|
@ -25,3 +25,4 @@ H2O-qmmm-gauss-16.inp 7
|
|||
H2O-qmmm-gauss-17.inp 7
|
||||
H2O-qmmm-gauss-18.inp 7
|
||||
H2O-qmmm-gauss-19.inp 7
|
||||
H2O-qmmm-none-force-mixing-1.inp 7
|
||||
|
|
|
|||
20
tests/QMMM/QS/sample_pdb/5H2O.pdb
Normal file
20
tests/QMMM/QS/sample_pdb/5H2O.pdb
Normal file
|
|
@ -0,0 +1,20 @@
|
|||
ATOM 1 O WAT 1 4.000 0.000 0.000 1.00 0.00 O
|
||||
ATOM 2 H1 WAT 1 4.000 0.000 0.947 1.00 0.00 H
|
||||
ATOM 3 H2 WAT 1 4.000 0.896 -0.317 1.00 0.00 H
|
||||
TER
|
||||
ATOM 4 O WAT 2 6.000 0.000 0.000 1.00 0.00 O
|
||||
ATOM 5 H1 WAT 2 6.000 0.000 0.947 1.00 0.00 H
|
||||
ATOM 6 H2 WAT 2 6.000 0.896 -0.317 1.00 0.00 H
|
||||
TER
|
||||
ATOM 7 O WAT 3 8.000 0.000 0.000 1.00 0.00 O
|
||||
ATOM 8 H1 WAT 3 8.000 0.000 0.947 1.00 0.00 H
|
||||
ATOM 9 H2 WAT 3 8.000 0.896 -0.317 1.00 0.00 H
|
||||
TER
|
||||
ATOM 10 O WAT 4 10.000 0.000 0.000 1.00 0.00 O
|
||||
ATOM 11 H1 WAT 4 10.000 0.000 0.947 1.00 0.00 H
|
||||
ATOM 12 H2 WAT 4 10.000 0.896 -0.317 1.00 0.00 H
|
||||
TER
|
||||
ATOM 13 O WAT 5 12.000 0.000 0.000 1.00 0.00 O
|
||||
ATOM 14 H1 WAT 5 12.000 0.000 0.947 1.00 0.00 H
|
||||
ATOM 15 H2 WAT 5 12.000 0.896 -0.317 1.00 0.00 H
|
||||
END
|
||||
76
tests/QMMM/QS/sample_pot/5H2O.pot
Normal file
76
tests/QMMM/QS/sample_pot/5H2O.pot
Normal file
|
|
@ -0,0 +1,76 @@
|
|||
*>>>>>>> AMBER FF Converted into CHARMM FF style <<<<<<<
|
||||
*>>>>>>> Generated on :: 20041222 165204.289 +0100 by :: teo <<<<<<<
|
||||
*>>>>>>> Teddy <<<<<<<
|
||||
*>>>>>>> Leap Title :: <<<<<<<
|
||||
*>>>>>>> Send all comments related to the FFs conversion to <<<<<<<
|
||||
*>>>>>>> teodoro.laino@gmail.com <<<<<<<
|
||||
|
||||
BONDS
|
||||
!
|
||||
!V(bond) = Kb(b - b0)**2
|
||||
!
|
||||
!Kb: kcal/mole/A**2
|
||||
!b0: A
|
||||
!
|
||||
!atom type Kb b0
|
||||
!
|
||||
HW OW 553.000000000 0.957200000
|
||||
HW HW 553.000000000 1.513600000
|
||||
|
||||
ANGLES
|
||||
!
|
||||
!V(angle) = Ktheta(Theta - Theta0)**2
|
||||
!
|
||||
!V(Urey-Bradley) = Kub(S - S0)**2
|
||||
!
|
||||
!Ktheta: kcal/mole/rad**2
|
||||
!Theta0: degrees
|
||||
!Kub: kcal/mole/A**2 (Urey-Bradley)
|
||||
!S0: A
|
||||
!
|
||||
!atom types Ktheta Theta0 Kub S0
|
||||
!
|
||||
|
||||
DIHEDRALS
|
||||
!
|
||||
!V(dihedral) = Kchi(1 + cos(n(chi) - delta))
|
||||
!
|
||||
!Kchi: kcal/mole
|
||||
!n: multiplicity
|
||||
!delta: degrees
|
||||
!
|
||||
!atom types Kchi n delta
|
||||
!
|
||||
|
||||
IMPROPER
|
||||
!
|
||||
!V(improper) = Kpsi(psi - psi0)**2
|
||||
!
|
||||
!Kpsi: kcal/mole/rad**2
|
||||
!psi0: degrees
|
||||
!note that the second column of numbers (0) is ignored
|
||||
!
|
||||
!atom types Kpsi psi0
|
||||
!
|
||||
|
||||
NONBONDED
|
||||
!
|
||||
!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
|
||||
!
|
||||
!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
|
||||
!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
|
||||
!
|
||||
!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
|
||||
!
|
||||
OW 0.000000000 0.152000000 1.768300000
|
||||
HW 0.000000000 0.000000000 0.000000000
|
||||
|
||||
END
|
||||
|
||||
!
|
||||
! This Section can be cutted & pasted into the Fist input file..
|
||||
!
|
||||
CHARGES
|
||||
OW -0.834000000
|
||||
HW 0.417000000
|
||||
END CHARGES
|
||||
45
tests/QMMM/QS/sample_psf/5H2O.psf
Normal file
45
tests/QMMM/QS/sample_psf/5H2O.psf
Normal file
|
|
@ -0,0 +1,45 @@
|
|||
PSF
|
||||
|
||||
1 !NTITLE
|
||||
Conversion from AMBER PARMTOP ::
|
||||
|
||||
15 !NATOM
|
||||
1 MOL1 1 WAT O OW -0.83400 16.00000 0
|
||||
2 MOL1 1 WAT H1 HW 0.41700 1.00800 0
|
||||
3 MOL1 1 WAT H2 HW 0.41700 1.00800 0
|
||||
4 MOL2 2 WAT O OW -0.83400 16.00000 0
|
||||
5 MOL2 2 WAT H1 HW 0.41700 1.00800 0
|
||||
6 MOL2 2 WAT H2 HW 0.41700 1.00800 0
|
||||
7 MOL3 3 WAT O OW -0.83400 16.00000 0
|
||||
8 MOL3 3 WAT H1 HW 0.41700 1.00800 0
|
||||
9 MOL3 3 WAT H2 HW 0.41700 1.00800 0
|
||||
10 MOL4 4 WAT O OW -0.83400 16.00000 0
|
||||
11 MOL4 4 WAT H1 HW 0.41700 1.00800 0
|
||||
12 MOL4 4 WAT H2 HW 0.41700 1.00800 0
|
||||
13 MOL5 5 WAT O OW -0.83400 16.00000 0
|
||||
14 MOL5 5 WAT H1 HW 0.41700 1.00800 0
|
||||
15 MOL5 5 WAT H2 HW 0.41700 1.00800 0
|
||||
|
||||
15 !NBOND
|
||||
3 1 2 1 2 3 6 4
|
||||
5 4 5 6 9 7 8 7
|
||||
8 9 10 11 10 12 11 12
|
||||
13 14 13 15 14 15
|
||||
|
||||
0 !NTHETA
|
||||
|
||||
|
||||
0 !NPHI
|
||||
|
||||
|
||||
0 !NIMPHI
|
||||
|
||||
|
||||
0 !NDON
|
||||
|
||||
0 !NACC
|
||||
|
||||
0 !NNB
|
||||
|
||||
0 !NGRP
|
||||
|
||||
136
tests/QMMM/QS/sample_xyz/5H2O_adaptive.xyz
Normal file
136
tests/QMMM/QS/sample_xyz/5H2O_adaptive.xyz
Normal file
|
|
@ -0,0 +1,136 @@
|
|||
15
|
||||
all intermolecular distances r = 2.0 <= r_inner = 2.5, should be 6/12
|
||||
O 5.000 0.000 0.000
|
||||
H 5.000 0.000 0.947
|
||||
H 5.000 0.896 -0.317
|
||||
O 7.000 0.000 0.000
|
||||
H 7.000 0.000 0.947
|
||||
H 7.000 0.896 -0.317
|
||||
O 9.000 0.000 0.000
|
||||
H 9.000 0.000 0.947
|
||||
H 9.000 0.896 -0.317
|
||||
O 11.000 0.000 0.000
|
||||
H 11.000 0.000 0.947
|
||||
H 11.000 0.896 -0.317
|
||||
O 13.000 0.000 0.000
|
||||
H 13.000 0.000 0.947
|
||||
H 13.000 0.896 -0.317
|
||||
15
|
||||
move molecule #4 to r_inner = 2.5 < r = 3.0 <= r_outer = 3.5, still 6/12
|
||||
O 5.100 0.000 0.000
|
||||
H 5.100 0.000 0.947
|
||||
H 5.100 0.896 -0.317
|
||||
O 7.100 0.000 0.000
|
||||
H 7.100 0.000 0.947
|
||||
H 7.100 0.896 -0.317
|
||||
O 9.100 0.000 0.000
|
||||
H 9.100 0.000 0.947
|
||||
H 9.100 0.896 -0.317
|
||||
O 12.100 0.000 0.000
|
||||
H 12.100 0.000 0.947
|
||||
H 12.100 0.896 -0.317
|
||||
O 14.100 0.000 0.000
|
||||
H 14.100 0.000 0.947
|
||||
H 14.100 0.896 -0.317
|
||||
15
|
||||
move 1 atom in molecule #4 to r = 3.55 > r_outer = 3.5, rest are at r = 3.45 <= r_outer, still 6/12
|
||||
O 5.100 0.000 0.000
|
||||
H 5.100 0.000 0.947
|
||||
H 5.100 0.896 -0.317
|
||||
O 7.100 0.000 0.000
|
||||
H 7.100 0.000 0.947
|
||||
H 7.100 0.896 -0.317
|
||||
O 9.100 0.000 0.000
|
||||
H 9.100 0.000 0.947
|
||||
H 9.100 0.896 -0.317
|
||||
O 12.650 0.000 0.000
|
||||
H 12.550 0.000 0.947
|
||||
H 12.550 0.896 -0.317
|
||||
O 14.650 0.000 0.000
|
||||
H 14.650 0.000 0.947
|
||||
H 14.650 0.896 -0.317
|
||||
15
|
||||
move rest of molecule #4 to r = 4.1 > r_outer = 3.5, now 6/9
|
||||
O 5.100 0.000 0.000
|
||||
H 5.100 0.000 0.947
|
||||
H 5.100 0.896 -0.317
|
||||
O 7.100 0.000 0.000
|
||||
H 7.100 0.000 0.947
|
||||
H 7.100 0.896 -0.317
|
||||
O 9.100 0.000 0.000
|
||||
H 9.100 0.000 0.947
|
||||
H 9.100 0.896 -0.317
|
||||
O 13.200 0.000 0.000
|
||||
H 13.200 0.000 0.947
|
||||
H 13.200 0.896 -0.317
|
||||
O 15.200 0.000 0.000
|
||||
H 15.200 0.000 0.947
|
||||
H 15.200 0.896 -0.317
|
||||
15
|
||||
move one atom of molecule #4 to r_inner = 2.5 < r = 3.45 < r_outer = 3.5, rest at r = 3.55 <= r_outer = 3.5 still 6/9
|
||||
O 5.100 0.000 0.000
|
||||
H 5.100 0.000 0.947
|
||||
H 5.100 0.896 -0.317
|
||||
O 7.100 0.000 0.000
|
||||
H 7.100 0.000 0.947
|
||||
H 7.100 0.896 -0.317
|
||||
O 9.100 0.000 0.000
|
||||
H 9.100 0.000 0.947
|
||||
H 9.100 0.896 -0.317
|
||||
O 12.550 0.000 0.000
|
||||
H 12.650 0.000 0.947
|
||||
H 12.650 0.896 -0.317
|
||||
O 14.650 0.000 0.000
|
||||
H 14.650 0.000 0.947
|
||||
H 14.650 0.896 -0.317
|
||||
15
|
||||
move rest of molecule #4 to r_inner = 2.5 < r = 3.1 < r_outer = 3.5, still 6/9
|
||||
O 5.100 0.000 0.000
|
||||
H 5.100 0.000 0.947
|
||||
H 5.100 0.896 -0.317
|
||||
O 7.100 0.000 0.000
|
||||
H 7.100 0.000 0.947
|
||||
H 7.100 0.896 -0.317
|
||||
O 9.100 0.000 0.000
|
||||
H 9.100 0.000 0.947
|
||||
H 9.100 0.896 -0.317
|
||||
O 12.200 0.000 0.000
|
||||
H 12.200 0.000 0.947
|
||||
H 12.200 0.896 -0.317
|
||||
O 14.200 0.000 0.000
|
||||
H 14.200 0.000 0.947
|
||||
H 14.200 0.896 -0.317
|
||||
15
|
||||
move one atom of molecule #4 to r = 2.45 < r_inner = 2.5, rest at r_inner = 2.5 < r = 2.55 <= r_outer = 3.5, back to 6/12
|
||||
O 5.100 0.000 0.000
|
||||
H 5.100 0.000 0.947
|
||||
H 5.100 0.896 -0.317
|
||||
O 7.100 0.000 0.000
|
||||
H 7.100 0.000 0.947
|
||||
H 7.100 0.896 -0.317
|
||||
O 9.100 0.000 0.000
|
||||
H 9.100 0.000 0.947
|
||||
H 9.100 0.896 -0.317
|
||||
O 11.550 0.000 0.000
|
||||
H 11.650 0.000 0.947
|
||||
H 11.650 0.896 -0.317
|
||||
O 13.650 0.000 0.000
|
||||
H 13.650 0.000 0.947
|
||||
H 13.650 0.896 -0.317
|
||||
15
|
||||
move rest of molecule #4 to r = 2.1 < r_inner = 2.5, still 6/12
|
||||
O 5.100 0.000 0.000
|
||||
H 5.100 0.000 0.947
|
||||
H 5.100 0.896 -0.317
|
||||
O 7.100 0.000 0.000
|
||||
H 7.100 0.000 0.947
|
||||
H 7.100 0.896 -0.317
|
||||
O 9.100 0.000 0.000
|
||||
H 9.100 0.000 0.947
|
||||
H 9.100 0.896 -0.317
|
||||
O 11.200 0.000 0.000
|
||||
H 11.200 0.000 0.947
|
||||
H 11.200 0.896 -0.317
|
||||
O 13.200 0.000 0.000
|
||||
H 13.200 0.000 0.947
|
||||
H 13.200 0.896 -0.317
|
||||
Loading…
Add table
Add a link
Reference in a new issue