From d4189efaa9ec1f12aca39a8bf713266fbc3aad9b Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Wed, 3 Mar 2021 17:07:25 -0500
Subject: [PATCH] sPuReMD: rework QMMM interface to use string representations
of atom element types.
---
sPuReMD/src/spuremd.c | 79 +++++++++++--------
sPuReMD/src/spuremd.h | 8 +-
.../qmmm_amber/fortran_stub/driver.F03 | 41 ++--------
.../qm2_extern_reaxff_puremd_module.F03 | 32 ++++----
4 files changed, 72 insertions(+), 88 deletions(-)
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index a39aa774..f4ab6109 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -37,6 +37,8 @@
#include "tool_box.h"
#include "vector.h"
+#include <ctype.h>
+
/* Handles additional entire geometry calculations after
* perturbing atom positions during a simulation step
@@ -954,10 +956,10 @@ int set_control_parameter( const void * const handle, const char * const keyword
* for the first simulation
*
* qm_num_atoms: num. atoms in the QM region
- * qm_types: element types for QM atoms
+ * qm_symbols: element types for QM atoms
* qm_pos: coordinates of QM atom positions (consecutively arranged), in Angstroms
* mm_num_atoms: num. atoms in the MM region
- * mm_types: element types for MM atoms
+ * mm_symbols: element types for MM atoms
* mm_pos_q: coordinates and charges of MM atom positions (consecutively arranged), in Angstroms / Coulombs
* sim_box_info: simulation box information, where the entries are
* - box length per dimension (3 entries)
@@ -965,13 +967,13 @@ int set_control_parameter( const void * const handle, const char * const keyword
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*/
-void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
- const double * const qm_pos, int mm_num_atoms, const int * const mm_types,
+void * setup_qmmm( int qm_num_atoms, const char * const qm_symbols,
+ const double * const qm_pos, int mm_num_atoms, const char * const mm_symbols,
const double * const mm_pos_q, const double * const sim_box_info,
const char * const ffield_file, const char * const control_file )
{
int i;
-// char atom_name[9];
+ char element[3];
rvec x;
spuremd_handle *spmd_handle;
@@ -1000,6 +1002,8 @@ void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
sim_box_info[3], sim_box_info[4], sim_box_info[5],
&spmd_handle->system->box );
+ element[2] = '\0';
+
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
x[0] = qm_pos[3 * i];
@@ -1009,12 +1013,14 @@ void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
Fit_to_Periodic_Box( &spmd_handle->system->box, x );
spmd_handle->workspace->orig_id[i] = i + 1;
-// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
-// element, sizeof(element) );
- spmd_handle->system->atoms[i].type = qm_types[i];
-// strncpy( spmd_handle->system->atoms[i].name, atom_name,
-// sizeof(spmd_handle->system->atoms[i].name) - 1 );
-// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ element[0] = toupper( qm_symbols[2 * i] );
+ element[1] = toupper( qm_symbols[2 * i + 1] );
+ Trim_Spaces( element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = Get_Atom_Type( &spmd_handle->system->reax_param,
+ element, sizeof(element) );
+ strncpy( spmd_handle->system->atoms[i].name, element,
+ sizeof(spmd_handle->system->atoms[i].name) - 1 );
+ spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
rvec_Copy( spmd_handle->system->atoms[i].x, x );
rvec_MakeZero( spmd_handle->system->atoms[i].v );
rvec_MakeZero( spmd_handle->system->atoms[i].f );
@@ -1033,12 +1039,14 @@ void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
Fit_to_Periodic_Box( &spmd_handle->system->box, x );
spmd_handle->workspace->orig_id[i] = i + 1;
-// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
-// element, sizeof(element) );
- spmd_handle->system->atoms[i].type = mm_types[i - spmd_handle->system->N_qm];
-// strncpy( spmd_handle->system->atoms[i].name, atom_name,
-// sizeof(spmd_handle->system->atoms[i].name) - 1 );
-// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ element[0] = toupper( mm_symbols[2 * (i - spmd_handle->system->N_qm)] );
+ element[1] = toupper( mm_symbols[2 * (i - spmd_handle->system->N_qm) + 1] );
+ Trim_Spaces( element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = Get_Atom_Type( &spmd_handle->system->reax_param,
+ element, sizeof(element) );
+ strncpy( spmd_handle->system->atoms[i].name, element,
+ sizeof(spmd_handle->system->atoms[i].name) - 1 );
+ spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
rvec_Copy( spmd_handle->system->atoms[i].x, x );
rvec_MakeZero( spmd_handle->system->atoms[i].v );
rvec_MakeZero( spmd_handle->system->atoms[i].f );
@@ -1059,10 +1067,10 @@ void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
*
* handle: pointer to wrapper struct with top-level data structures
* qm_num_atoms: num. atoms in the QM region
- * qm_types: element types for QM atoms
+ * qm_symbols: element types for QM atoms
* qm_pos: coordinates of QM atom positions (consecutively arranged), in Angstroms
* mm_num_atoms: num. atoms in the MM region
- * mm_types: element types for MM atoms
+ * mm_symbols: element types for MM atoms
* mm_pos_q: coordinates and charges of MM atom positions (consecutively arranged), in Angstroms / Coulombs
* sim_box_info: simulation box information, where the entries are
* - box length per dimension (3 entries)
@@ -1073,12 +1081,13 @@ void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int reset_qmmm( const void * const handle, int qm_num_atoms,
- const int * const qm_types, const double * const qm_pos,
- int mm_num_atoms, const int * const mm_types,
+ const char * const qm_symbols, const double * const qm_pos,
+ int mm_num_atoms, const char * const mm_symbols,
const double * const mm_pos_q, const double * const sim_box_info,
const char * const ffield_file, const char * const control_file )
{
int i, ret;
+ char element[3];
rvec x;
spuremd_handle *spmd_handle;
@@ -1127,12 +1136,14 @@ int reset_qmmm( const void * const handle, int qm_num_atoms,
Fit_to_Periodic_Box( &spmd_handle->system->box, x );
spmd_handle->workspace->orig_id[i] = i + 1;
-// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
-// element, sizeof(element) );
- spmd_handle->system->atoms[i].type = qm_types[i];
-// strncpy( spmd_handle->system->atoms[i].name, atom_name,
-// sizeof(spmd_handle->system->atoms[i].name) - 1 );
-// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ element[0] = toupper( qm_symbols[2 * i] );
+ element[1] = toupper( qm_symbols[2 * i + 1] );
+ Trim_Spaces( element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = Get_Atom_Type( &spmd_handle->system->reax_param,
+ element, sizeof(element) );
+ strncpy( spmd_handle->system->atoms[i].name, element,
+ sizeof(spmd_handle->system->atoms[i].name) - 1 );
+ spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
rvec_Copy( spmd_handle->system->atoms[i].x, x );
rvec_MakeZero( spmd_handle->system->atoms[i].v );
rvec_MakeZero( spmd_handle->system->atoms[i].f );
@@ -1150,12 +1161,14 @@ int reset_qmmm( const void * const handle, int qm_num_atoms,
Fit_to_Periodic_Box( &spmd_handle->system->box, x );
spmd_handle->workspace->orig_id[i] = i + 1;
-// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
-// element, sizeof(element) );
- spmd_handle->system->atoms[i].type = mm_types[i - spmd_handle->system->N_qm];
-// strncpy( spmd_handle->system->atoms[i].name, atom_name,
-// sizeof(spmd_handle->system->atoms[i].name) - 1 );
-// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ element[0] = toupper( mm_symbols[2 * (i - spmd_handle->system->N_qm)] );
+ element[1] = toupper( mm_symbols[2 * (i - spmd_handle->system->N_qm) + 1] );
+ Trim_Spaces( element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = Get_Atom_Type( &spmd_handle->system->reax_param,
+ element, sizeof(element) );
+ strncpy( spmd_handle->system->atoms[i].name, element,
+ sizeof(spmd_handle->system->atoms[i].name) - 1 );
+ spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
rvec_Copy( spmd_handle->system->atoms[i].x, x );
rvec_MakeZero( spmd_handle->system->atoms[i].v );
rvec_MakeZero( spmd_handle->system->atoms[i].f );
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index cae887f4..22f534e5 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -73,13 +73,13 @@ int set_control_parameter( const void * const, const char * const,
const char ** const );
#if defined(QMMM)
-void * setup_qmmm( int, const int * const,
- const double * const, int, const int * const,
+void * setup_qmmm( int, const char * const,
+ const double * const, int, const char * const,
const double * const, const double * const,
const char * const, const char * const );
-int reset_qmmm( const void * const, int, const int * const,
- const double * const, int, const int * const,
+int reset_qmmm( const void * const, int, const char * const,
+ const double * const, int, const char * const,
const double * const, const double * const,
const char * const, const char * const );
diff --git a/test/library/qmmm_amber/fortran_stub/driver.F03 b/test/library/qmmm_amber/fortran_stub/driver.F03
index 7532f50e..68a287cf 100644
--- a/test/library/qmmm_amber/fortran_stub/driver.F03
+++ b/test/library/qmmm_amber/fortran_stub/driver.F03
@@ -14,7 +14,7 @@ contains
integer (c_int), intent(in) :: num_qm_atoms, num_mm_atoms
real (c_double), intent(out) :: sim_box_info(6)
real (c_double), intent(out) :: qm_pos(3*num_qm_atoms), qm_q(num_qm_atoms), mm_pos_q(4*num_mm_atoms)
- character(len=5), intent(out) :: qm_types(num_qm_atoms), mm_types(num_mm_atoms)
+ character(len=2), intent(out) :: qm_types(num_qm_atoms), mm_types(num_mm_atoms)
character(len=200) :: line
integer :: bgf_version, num_descrp, num_remark, num_format, line_num
@@ -59,7 +59,7 @@ contains
if (line(1:6) == 'HETATM') then
if (bgf_version < 400) then
if (nqm < num_qm_atoms) then
- read (line, '(30x,3f10.5,1x,a5,6x,f10.5)', end=40, err=40) &
+ read (line, '(30x,3f10.5,4x,a2,6x,f10.5)', end=40, err=40) &
qm_pos(3*nqm+1), qm_pos(3*nqm+2), qm_pos(3*nqm+3), &
qm_types(nqm+1), qm_q(nqm+1)
@@ -67,7 +67,7 @@ contains
else
! take the MM charges from BGF file;
! Amber would set instead in real QM/MM simulation
- read (line, '(30x,3f10.5,1x,a5,6x,f10.5)', end=40, err=40) &
+ read (line, '(30x,3f10.5,4x,a2,6x,f10.5)', end=40, err=40) &
mm_pos_q(4*nmm+1), mm_pos_q(4*nmm+2), mm_pos_q(4*nmm+3), &
mm_types(nmm+1), mm_pos_q(4*nmm+4)
@@ -117,8 +117,7 @@ program driver
integer (c_int) :: num_qm_atoms, num_mm_atoms
real (c_double) :: sim_box_info(6), e_total
real (c_double), dimension(:), allocatable :: qm_pos, qm_q, qm_f, mm_pos_q, mm_f
- character(len=5), dimension(:), allocatable :: qm_types, mm_types
- integer (c_int), dimension(:), allocatable :: qm_types_map, mm_types_map
+ character(len=2), dimension(:), allocatable :: qm_types, mm_types
num_args = command_argument_count( )
if (num_args /= 3) then
@@ -138,12 +137,10 @@ program driver
read (args(3), *) num_mm_atoms
allocate( qm_types(num_qm_atoms) )
- allocate( qm_types_map(num_qm_atoms) )
allocate( qm_pos(3*num_qm_atoms) )
allocate( qm_q(num_qm_atoms) )
allocate( qm_f(3*num_qm_atoms) )
allocate( mm_types(num_mm_atoms) )
- allocate( mm_types_map(num_mm_atoms) )
allocate( mm_pos_q(4*num_mm_atoms) )
allocate( mm_f(3*num_mm_atoms) )
@@ -158,32 +155,8 @@ program driver
stop
endif
- do ix = 1, num_qm_atoms
- if (trim(adjustl(qm_types(ix))) == 'C') then
- qm_types_map(ix) = 0
- else if (trim(adjustl(qm_types(ix))) == 'H') then
- qm_types_map(ix) = 1
- else if (trim(adjustl(qm_types(ix))) == 'O') then
- qm_types_map(ix) = 2
- else
- stop 'ERROR: unrecogized QM atom type'
- endif
- end do
-
- do ix = 1, num_mm_atoms
- if (trim(adjustl(mm_types(ix))) == 'C') then
- mm_types_map(ix) = 0
- else if (trim(adjustl(mm_types(ix))) == 'H') then
- mm_types_map(ix) = 1
- else if (trim(adjustl(mm_types(ix))) == 'O') then
- mm_types_map(ix) = 2
- else
- stop 'ERROR: unrecogized MM atom type'
- endif
- end do
-
- call get_reaxff_puremd_forces( num_qm_atoms, qm_pos, qm_types_map, &
- qm_q, num_mm_atoms, mm_pos_q, mm_types_map, e_total, qm_f, &
+ call get_reaxff_puremd_forces( num_qm_atoms, qm_pos, qm_types, &
+ qm_q, num_mm_atoms, mm_pos_q, mm_types, e_total, qm_f, &
mm_f, qmcharge )
write (stdout, fmt=75) 'Total energy:', e_total
@@ -202,12 +175,10 @@ program driver
deallocate(args)
deallocate(qm_types)
- deallocate(qm_types_map)
deallocate(qm_pos)
deallocate(qm_q)
deallocate(qm_f)
deallocate(mm_types)
- deallocate(mm_types_map)
deallocate(mm_pos_q)
deallocate(mm_f)
diff --git a/test/library/qmmm_amber/fortran_stub/qm2_extern_reaxff_puremd_module.F03 b/test/library/qmmm_amber/fortran_stub/qm2_extern_reaxff_puremd_module.F03
index 772c5354..8f7469e2 100644
--- a/test/library/qmmm_amber/fortran_stub/qm2_extern_reaxff_puremd_module.F03
+++ b/test/library/qmmm_amber/fortran_stub/qm2_extern_reaxff_puremd_module.F03
@@ -5,17 +5,17 @@ module qm2_extern_reaxff_puremd_module
interface
type(c_ptr) function setup_qmmm &
- (num_qm_atoms, qm_types, qm_pos, &
- num_mm_atoms, mm_types, mm_pos_q, sim_box_info, &
+ (num_qm_atoms, qm_symbols, qm_pos, &
+ num_mm_atoms, mm_symbols, mm_pos_q, sim_box_info, &
ffield_filename, control_filename) &
bind(C, name='setup_qmmm')
use, intrinsic :: iso_c_binding
implicit none
integer (c_int), value :: num_qm_atoms
- type(c_ptr), value :: qm_types
+ type(c_ptr), value :: qm_symbols
type(c_ptr), value :: qm_pos
integer (c_int), value :: num_mm_atoms
- type(c_ptr), value :: mm_types
+ type(c_ptr), value :: mm_symbols
type(c_ptr), value :: mm_pos_q
type(c_ptr), value :: sim_box_info
type(c_ptr), value :: ffield_filename
@@ -23,18 +23,18 @@ module qm2_extern_reaxff_puremd_module
end function setup_qmmm
integer (c_int) function reset_qmmm &
- (handle, num_qm_atoms, qm_types, qm_pos, &
- num_mm_atoms, mm_types, mm_pos_q, sim_box_info, &
+ (handle, num_qm_atoms, qm_symbols, qm_pos, &
+ num_mm_atoms, mm_symbols, mm_pos_q, sim_box_info, &
ffield_filename, control_filename) &
bind(C, name='reset_qmmm')
use, intrinsic :: iso_c_binding
implicit none
type(c_ptr), value :: handle
integer (c_int), value :: num_qm_atoms
- type(c_ptr), value :: qm_types
+ type(c_ptr), value :: qm_symbols
type(c_ptr), value :: qm_pos
integer (c_int), value :: num_mm_atoms
- type(c_ptr), value :: mm_types
+ type(c_ptr), value :: mm_symbols
type(c_ptr), value :: mm_pos_q
type(c_ptr), value :: sim_box_info
type(c_ptr), value :: ffield_filename
@@ -119,19 +119,19 @@ module qm2_extern_reaxff_puremd_module
type(c_ptr), save, private :: handle = c_null_ptr
contains
- subroutine get_reaxff_puremd_forces( num_qm_atoms, qm_pos, qm_types, &
- qm_q, num_mm_atoms, mm_pos_q, mm_types, e_total, &
+ subroutine get_reaxff_puremd_forces( num_qm_atoms, qm_pos, qm_symbols, &
+ qm_q, num_mm_atoms, mm_pos_q, mm_symbols, e_total, &
qm_f, mm_f, qmcharge )
use, intrinsic :: iso_c_binding
implicit none
integer (c_int), intent(in) :: num_qm_atoms ! number of QM atoms
real (c_double), target, intent(in) :: qm_pos(3,num_qm_atoms) ! QM atom coordinates
- integer (c_int), target, intent(in) :: qm_types(num_qm_atoms) ! QM atom types
+ character(len=2), dimension(num_qm_atoms), target, intent(in) :: qm_symbols ! QM atom types
real (c_double), target, intent(inout) :: qm_q(num_qm_atoms) ! QM atom charges (nuclear charge in au)
integer (c_int), intent(in) :: num_mm_atoms ! number of MM atoms
real (c_double), target, intent(in) :: mm_pos_q(4,num_mm_atoms) ! MM atom coordinates and charges (nuclear charge in au)
- integer (c_int), target, intent(in) :: mm_types(num_mm_atoms) ! MM atom types
+ character(len=2), dimension(num_qm_atoms), target, intent(in) :: mm_symbols ! MM atom types
real (c_double), target, intent(out) :: e_total ! SCF energy (kcal/mol)
real (c_double), target, intent(out) :: qm_f(3,num_qm_atoms) ! SCF QM force (AMU * Angstroms / ps^2)
real (c_double), target, intent(out) :: mm_f(3,num_mm_atoms) ! SCF MM force (AMU * Angstroms / ps^2)
@@ -159,8 +159,8 @@ contains
if ( first_call ) then
first_call = .false.
- handle = setup_qmmm( num_qm_atoms, c_loc(qm_types), &
- c_loc(qm_pos), num_mm_atoms, c_loc(mm_types), &
+ handle = setup_qmmm( num_qm_atoms, c_loc(qm_symbols), &
+ c_loc(qm_pos), num_mm_atoms, c_loc(mm_symbols), &
c_loc(mm_pos_q), c_loc(sim_box_info), c_loc(ffield_filename), c_null_ptr )
! NVE ensemble
@@ -297,8 +297,8 @@ contains
c_loc(e_total), c_null_ptr, c_null_ptr, c_null_ptr )
if ( ret /= 0_c_int ) stop "ERROR: get_reaxff_puremd_forces::get_system_info"
else
- ret = reset_qmmm( handle, num_qm_atoms, c_loc(qm_types), &
- c_loc(qm_pos), num_mm_atoms, c_loc(mm_types), &
+ ret = reset_qmmm( handle, num_qm_atoms, c_loc(qm_symbols), &
+ c_loc(qm_pos), num_mm_atoms, c_loc(mm_symbols), &
c_loc(mm_pos_q), c_loc(sim_box_info), c_null_ptr, c_null_ptr )
if ( ret /= 0_c_int ) stop "ERROR: get_reaxff_puremd_forces::reset_qmmm"
--
GitLab