-
Kurt A. O'Hearn authored
sPuReMD: rework Fortran F90 interfaces for Amber specifics and declare this interface for such usage.
Kurt A. O'Hearn authoredsPuReMD: rework Fortran F90 interfaces for Amber specifics and declare this interface for such usage.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
spuremd.c 54.68 KiB
/*----------------------------------------------------------------------
SerialReax - Reax Force Field Simulator
Copyright (2010) Purdue University
Hasan Metin Aktulga, haktulga@cs.purdue.edu
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details:
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "spuremd.h"
#include "allocate.h"
#include "analyze.h"
#include "box.h"
#include "control.h"
#include "ffield.h"
#include "forces.h"
#include "init_md.h"
#include "io_tools.h"
#include "neighbors.h"
#include "geo_tools.h"
#include "reset_tools.h"
#include "restart.h"
#include "system_props.h"
#include "tool_box.h"
#include "vector.h"
/* Handles additional entire geometry calculations after
* perturbing atom positions during a simulation step
*/
static void Post_Evolve( reax_system * const system, control_params * const control,
simulation_data * const data, static_storage * const workspace,
reax_list ** const lists, output_controls * const out_control )
{
int i;
rvec diff, cross;
/* remove rotational and translational velocity of the center of mass */
if ( control->ensemble != NVE && control->remove_CoM_vel > 0
&& data->step % control->remove_CoM_vel == 0 )
{
/* compute velocity of the center of mass */
Compute_Center_of_Mass( system, data );
for ( i = 0; i < system->N; i++ )
{
/* remove translational */
rvec_ScaledAdd( system->atoms[i].v, -1.0, data->vcm );
/* remove rotational */
rvec_ScaledSum( diff, 1.0, system->atoms[i].x, -1.0, data->xcm );
rvec_Cross( cross, data->avcm, diff );
rvec_ScaledAdd( system->atoms[i].v, -1.0, cross );
}
}
if ( control->ensemble == NVE ) {
Compute_Kinetic_Energy( system, data );
}
if ( (out_control->log_update_freq > 0
&& data->step % out_control->log_update_freq == 0)
|| (out_control->write_steps > 0
&& data->step % out_control->write_steps == 0) )
{
Compute_Total_Energy( data );
}
if ( control->compute_pressure == TRUE && control->ensemble != sNPT
&& control->ensemble != iNPT && control->ensemble != aNPT )
{
Compute_Pressure_Isotropic( system, control, data, out_control );
}
}
/* Parse input files
*
* geo_file: file containing geometry info of the structure to simulate
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*/
static void Read_Input_Files( const char * const geo_file,
const char * const ffield_file, const char * const control_file,
reax_system * const system, control_params * const control,
simulation_data * const data, static_storage * const workspace,
output_controls * const out_control )
{
if ( ffield_file != NULL )
{
Read_Force_Field( ffield_file, system, &system->reax_param );
}
Set_Control_Defaults( system, control, out_control );
if ( control_file != NULL )
{
Read_Control_File( control_file, system, control, out_control );
}
Set_Control_Derived_Values( system, control );
if ( geo_file != NULL )
{
if ( control->geo_format == CUSTOM )
{
Read_Geo( geo_file, system, control, data, workspace );
}
else if ( control->geo_format == PDB )
{
Read_PDB( geo_file, system, control, data, workspace );
}
else if ( control->geo_format == BGF )
{
Read_BGF( geo_file, system, control, data, workspace );
}
else if ( control->geo_format == ASCII_RESTART )
{
Read_ASCII_Restart( geo_file, system, control, data, workspace );
control->restart = TRUE;
}
else if ( control->geo_format == BINARY_RESTART )
{
Read_Binary_Restart( geo_file, system, control, data, workspace );
control->restart = TRUE;
} else
{
fprintf( stderr, "[ERROR] unknown geo file format. terminating!\n" );
exit( INVALID_GEO );
}
}
#if defined(DEBUG_FOCUS)
Print_Box( &system->box, stderr );
#endif
}
/* Allocate top-level data structures and parse input files
* for the first simulation
*
* geo_file: file containing geometry info of the structure to simulate
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*/
void * setup( const char * const geo_file, const char * const ffield_file,
const char * const control_file )
{
int i;
spuremd_handle *spmd_handle;
/* top-level allocation */
spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
"setup::spmd_handle" );
/* second-level allocations */
spmd_handle->system = smalloc( sizeof(reax_system),
"Setup::spmd_handle->system" );
spmd_handle->system->prealloc_allocated = FALSE;
spmd_handle->system->ffield_params_allocated = FALSE;
spmd_handle->system->g.allocated = FALSE;
spmd_handle->control = smalloc( sizeof(control_params),
"Setup::spmd_handle->control" );
spmd_handle->data = smalloc( sizeof(simulation_data),
"Setup::spmd_handle->data" );
spmd_handle->workspace = smalloc( sizeof(static_storage),
"Setup::spmd_handle->workspace" );
spmd_handle->workspace->H.allocated = FALSE;
spmd_handle->workspace->H_full.allocated = FALSE;
spmd_handle->workspace->H_sp.allocated = FALSE;
spmd_handle->workspace->H_p.allocated = FALSE;
spmd_handle->workspace->H_spar_patt.allocated = FALSE;
spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
spmd_handle->workspace->H_app_inv.allocated = FALSE;
spmd_handle->workspace->L.allocated = FALSE;
spmd_handle->workspace->U.allocated = FALSE;
spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
"Setup::spmd_handle->lists" );
for ( i = 0; i < LIST_N; ++i )
{
spmd_handle->lists[i] = smalloc( sizeof(reax_list),
"Setup::spmd_handle->lists[i]" );
spmd_handle->lists[i]->allocated = FALSE;
}
spmd_handle->out_control = smalloc( sizeof(output_controls),
"Setup::spmd_handle->out_control" );
spmd_handle->output_enabled = TRUE;
spmd_handle->realloc = TRUE;
spmd_handle->callback = NULL;
spmd_handle->data->sim_id = 0;
Read_Input_Files( geo_file, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
return (void *) spmd_handle;
}
/* Setup callback function to be run after each simulation step
*
* handle: pointer to wrapper struct with top-level data structures
* callback: function pointer to attach for callback
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int setup_callback( const void * const handle, const callback_function callback )
{
int ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL && callback != NULL )
{
spmd_handle = (spuremd_handle*) handle;
spmd_handle->callback = callback;
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Run the simulation according to the prescribed parameters
*
* handle: pointer to wrapper struct with top-level data structures
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int simulate( const void * const handle )
{
int steps, ret;
evolve_function Evolve;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
Initialize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists,
spmd_handle->out_control, &Evolve,
spmd_handle->output_enabled,
spmd_handle->realloc );
spmd_handle->realloc = FALSE;
/* compute f_0 */
//if( control.restart == FALSE ) {
Reset( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists );
Generate_Neighbor_Lists( spmd_handle->system, spmd_handle->control, spmd_handle->data, spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->realloc );
Compute_Kinetic_Energy( spmd_handle->system, spmd_handle->data );
if ( spmd_handle->control->compute_pressure == TRUE && spmd_handle->control->ensemble != sNPT
&& spmd_handle->control->ensemble != iNPT && spmd_handle->control->ensemble != aNPT )
{
Compute_Pressure_Isotropic( spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->out_control );
}
if ( spmd_handle->output_enabled == TRUE || spmd_handle->callback != NULL )
{
if ( ((spmd_handle->out_control->log_update_freq > 0
&& spmd_handle->data->step % spmd_handle->out_control->log_update_freq == 0)
|| (spmd_handle->out_control->write_steps > 0
&& spmd_handle->data->step % spmd_handle->out_control->write_steps == 0))
|| spmd_handle->callback != NULL )
{
Compute_Total_Energy( spmd_handle->data );
}
}
if ( spmd_handle->output_enabled == TRUE )
{
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
}
Check_Energy( spmd_handle->data );
if ( spmd_handle->output_enabled == TRUE )
{
if ( spmd_handle->out_control->write_steps > 0
&& spmd_handle->data->step % spmd_handle->out_control->write_steps == 0 )
{
Write_PDB( spmd_handle->system, spmd_handle->lists[BONDS], spmd_handle->data,
spmd_handle->control, spmd_handle->workspace, spmd_handle->out_control );
}
}
if ( spmd_handle->callback != NULL )
{
spmd_handle->callback( spmd_handle->system->N, spmd_handle->system->atoms,
spmd_handle->data );
}
//}
for ( ++spmd_handle->data->step; spmd_handle->data->step <= spmd_handle->control->nsteps; spmd_handle->data->step++ )
{
if ( spmd_handle->control->T_mode != 0 )
{
Temperature_Control( spmd_handle->control, spmd_handle->data,
spmd_handle->out_control );
}
Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
Post_Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
if ( spmd_handle->output_enabled == TRUE )
{
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control ); }
Check_Energy( spmd_handle->data );
if ( spmd_handle->output_enabled == TRUE )
{
steps = spmd_handle->data->step - spmd_handle->data->prev_steps;
Analysis( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
if ( spmd_handle->out_control->restart_freq > 0
&& steps % spmd_handle->out_control->restart_freq == 0 )
{
Write_Restart( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->out_control );
}
if ( spmd_handle->out_control->write_steps > 0
&& steps % spmd_handle->out_control->write_steps == 0 )
{
Write_PDB( spmd_handle->system, spmd_handle->lists[BONDS], spmd_handle->data,
spmd_handle->control, spmd_handle->workspace, spmd_handle->out_control );
}
}
if ( spmd_handle->callback != NULL )
{
spmd_handle->callback( spmd_handle->system->N, spmd_handle->system->atoms,
spmd_handle->data );
}
}
spmd_handle->data->timing.end = Get_Time( );
spmd_handle->data->timing.elapsed = Get_Timing_Info( spmd_handle->data->timing.start );
if ( spmd_handle->output_enabled == TRUE
&& spmd_handle->out_control->log_update_freq > 0 )
{
fprintf( spmd_handle->out_control->log, "total: %.2f secs\n", spmd_handle->data->timing.elapsed );
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Deallocate all data structures post-simulation
*
* handle: pointer to wrapper struct with top-level data structures
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int cleanup( const void * const handle )
{
int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->output_enabled, FALSE );
sfree( spmd_handle->out_control, "cleanup::spmd_handle->out_control" );
for ( i = 0; i < LIST_N; ++i )
{
sfree( spmd_handle->lists[i], "cleanup::spmd_handle->lists[i]" );
}
sfree( spmd_handle->lists, "cleanup::spmd_handle->lists" );
sfree( spmd_handle->workspace, "cleanup::spmd_handle->workspace" );
sfree( spmd_handle->data, "cleanup::spmd_handle->data" );
sfree( spmd_handle->control, "cleanup::spmd_handle->control" );
sfree( spmd_handle->system, "cleanup::spmd_handle->system" );
sfree( spmd_handle, "cleanup::spmd_handle" );
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Reset for the next simulation by parsing input files and triggering
* reallocation if more space is needed
*
* handle: pointer to wrapper struct with top-level data structures
* geo_file: file containing geometry info of the structure to simulate
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int reset( const void * const handle, const char * const geo_file,
const char * const ffield_file, const char * const control_file )
{
int ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
/* close files used in previous simulation */
if ( spmd_handle->output_enabled == TRUE )
{
Finalize_Out_Controls( spmd_handle->system, spmd_handle->control,
spmd_handle->workspace, spmd_handle->out_control );
}
spmd_handle->realloc = FALSE;
spmd_handle->data->sim_id++;
Read_Input_Files( geo_file, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
if ( spmd_handle->system->N > spmd_handle->system->N_max )
{
/* deallocate everything which needs more space
* (i.e., structures whose space is a function of the number of atoms),
* except for data structures allocated while parsing input files */
Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->output_enabled, TRUE );
spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
spmd_handle->realloc = TRUE;
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for atom positions
*
* handle: pointer to wrapper struct with top-level data structures
* pos: coordinates of atom positions, in Angstroms (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_positions( const void * const handle, double * const pos )
{
int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
for ( i = 0; i < spmd_handle->system->N; ++i )
{
pos[3 * i] = spmd_handle->system->atoms[i].x[0];
pos[3 * i + 1] = spmd_handle->system->atoms[i].x[1];
pos[3 * i + 2] = spmd_handle->system->atoms[i].x[2];
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for atom velocities
*
* handle: pointer to wrapper struct with top-level data structures
* vel: coordinates of atom velocities, in Angstroms / ps (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_velocities( const void * const handle, double * const vel )
{
int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
for ( i = 0; i < spmd_handle->system->N; ++i )
{
vel[3 * i] = spmd_handle->system->atoms[i].v[0];
vel[3 * i + 1] = spmd_handle->system->atoms[i].v[1];
vel[3 * i + 2] = spmd_handle->system->atoms[i].v[2];
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for atom forces
*
* handle: pointer to wrapper struct with top-level data structures
* f: coordinates of atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_forces( const void * const handle, double * const f )
{
int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
for ( i = 0; i < spmd_handle->system->N; ++i )
{
f[3 * i] = spmd_handle->system->atoms[i].f[0];
f[3 * i + 1] = spmd_handle->system->atoms[i].f[1];
f[3 * i + 2] = spmd_handle->system->atoms[i].f[2];
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for atom charges
*
* handle: pointer to wrapper struct with top-level data structures
* q: atom charges, in Coulombs (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_charges( const void * const handle, double * const q )
{
int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
for ( i = 0; i < spmd_handle->system->N; ++i )
{
q[i] = spmd_handle->system->atoms[i].q;
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for system energies
*
* handle: pointer to wrapper struct with top-level data structures
* e_pot: system potential energy, in kcal / mol (reference from caller)
* e_kin: system kinetic energy, in kcal / mol (reference from caller)
* e_tot: system total energy, in kcal / mol (reference from caller)
* t_scalar: temperature scalar, in K (reference from caller) * vol: volume of the simulation box, in Angstroms^3 (reference from caller)
* pres: average pressure, in K (reference from caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_system_info( const void * const handle, double * const e_pot,
double * const e_kin, double * const e_tot, double * const temp,
double * const vol, double * const pres )
{
int ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
if ( e_pot != NULL )
{
*e_pot = spmd_handle->data->E_Pot;
}
if ( e_kin != NULL )
{
*e_kin = spmd_handle->data->E_Kin;
}
if ( e_tot != NULL )
{
*e_tot = spmd_handle->data->E_Tot;
}
if ( temp != NULL )
{
*temp = spmd_handle->data->therm.T;
}
if ( vol != NULL )
{
*vol = spmd_handle->system->box.volume;
}
if ( pres != NULL )
{
*pres = (spmd_handle->control->P[0] + spmd_handle->control->P[1]
+ spmd_handle->control->P[2]) / 3.0;
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Setter for writing output to files
*
* handle: pointer to wrapper struct with top-level data structures
* enabled: TRUE enables writing output to files, FALSE otherwise
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int set_output_enabled( const void * const handle, const int enabled )
{
int ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
spmd_handle->output_enabled = enabled;
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Setter for simulation parameter values as defined in the input control file
*
* handle: pointer to wrapper struct with top-level data structures
* control_keyword: keyword from the control file to set the value for
* control_value: value to set
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int set_control_parameter( const void * const handle, const char * const keyword,
const char ** const values )
{
int ret, ret_;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
ret_ = Set_Control_Parameter( keyword, values, spmd_handle->control,
spmd_handle->out_control );
if ( ret_ == SUCCESS )
{
ret = SPUREMD_SUCCESS;
}
}
return ret;
}
#if defined(QMMM)
/* Allocate top-level data structures and parse input files
* for the first simulation
*
* qm_num_atoms: num. atoms in the QM region
* qm_types: 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_pos: coordinates of MM atom positions (consecutively arranged), in Angstroms
* mm_q: charge of MM atom, in Coulombs
* sim_box_info: simulation box information, where the entries are
* - box length per dimension (3 entries)
* - angles per dimension (3 entries)
* 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,
const double * const mm_pos, const double * const mm_q,
const double * const sim_box_info, const char * const ffield_file,
const char * const control_file )
{
int i;
// char atom_name[9];
rvec x;
spuremd_handle *spmd_handle;
/* top-level allocation */
spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
"setup::spmd_handle" );
/* second-level allocations */
spmd_handle->system = smalloc( sizeof(reax_system),
"Setup::spmd_handle->system" );
spmd_handle->system->prealloc_allocated = FALSE;
spmd_handle->system->ffield_params_allocated = FALSE;
spmd_handle->system->g.allocated = FALSE;
spmd_handle->control = smalloc( sizeof(control_params),
"Setup::spmd_handle->control" );
spmd_handle->data = smalloc( sizeof(simulation_data),
"Setup::spmd_handle->data" );
spmd_handle->workspace = smalloc( sizeof(static_storage),
"Setup::spmd_handle->workspace" );
spmd_handle->workspace->H.allocated = FALSE;
spmd_handle->workspace->H_full.allocated = FALSE;
spmd_handle->workspace->H_sp.allocated = FALSE;
spmd_handle->workspace->H_p.allocated = FALSE;
spmd_handle->workspace->H_spar_patt.allocated = FALSE;
spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
spmd_handle->workspace->H_app_inv.allocated = FALSE;
spmd_handle->workspace->L.allocated = FALSE;
spmd_handle->workspace->U.allocated = FALSE;
spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
"Setup::spmd_handle->lists" );
for ( i = 0; i < LIST_N; ++i )
{
spmd_handle->lists[i] = smalloc( sizeof(reax_list),
"Setup::spmd_handle->lists[i]" );
spmd_handle->lists[i]->allocated = FALSE;
}
spmd_handle->out_control = smalloc( sizeof(output_controls),
"Setup::spmd_handle->out_control" );
spmd_handle->output_enabled = FALSE;
spmd_handle->realloc = TRUE;
spmd_handle->callback = NULL;
spmd_handle->data->sim_id = 0;
spmd_handle->system->N_qm = qm_num_atoms;
spmd_handle->system->N_mm = mm_num_atoms;
spmd_handle->system->N = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
PreAllocate_Space( spmd_handle->system, spmd_handle->control,
spmd_handle->workspace, spmd_handle->system->N );
Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
sim_box_info[3], sim_box_info[4], sim_box_info[5],
&spmd_handle->system->box );
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
x[0] = qm_pos[3 * i];
x[1] = qm_pos[3 * i + 1];
x[2] = qm_pos[3 * i + 2];
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';
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 );
spmd_handle->system->atoms[i].q = 0.0;
spmd_handle->system->atoms[i].q_init = 0.0;
spmd_handle->system->atoms[i].qmmm_mask = TRUE;
}
for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
{
x[0] = mm_pos[3 * (i - spmd_handle->system->N_qm)];
x[1] = mm_pos[3 * (i - spmd_handle->system->N_qm) + 1];
x[2] = mm_pos[3 * (i - spmd_handle->system->N_qm) + 2];
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';
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 );
spmd_handle->system->atoms[i].q = mm_q[i - spmd_handle->system->N_qm];
spmd_handle->system->atoms[i].q_init = mm_q[i - spmd_handle->system->N_qm];
spmd_handle->system->atoms[i].qmmm_mask = FALSE;
}
Read_Input_Files( NULL, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
return (void *) spmd_handle;
}
/* Reset for the next simulation by parsing input files and triggering
* reallocation if more space is needed
*
* 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_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_pos: coordinates of MM atom positions (consecutively arranged), in Angstroms
* mm_q: charge of MM atom, in Coulombs
* sim_box_info: simulation box information, where the entries are
* - box length per dimension (3 entries)
* - angles per dimension (3 entries)
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*
* 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 double * const mm_pos, const double * const mm_q,
const double * const sim_box_info, const char * const ffield_file, const char * const control_file )
{
int i, ret;
rvec x;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
/* close files used in previous simulation */
if ( spmd_handle->output_enabled == TRUE )
{
Finalize_Out_Controls( spmd_handle->system, spmd_handle->control,
spmd_handle->workspace, spmd_handle->out_control );
}
spmd_handle->realloc = FALSE;
spmd_handle->data->sim_id++;
spmd_handle->system->N_qm = qm_num_atoms;
spmd_handle->system->N_mm = mm_num_atoms;
spmd_handle->system->N = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
PreAllocate_Space( spmd_handle->system, spmd_handle->control,
spmd_handle->workspace, spmd_handle->system->N );
Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
sim_box_info[3], sim_box_info[4], sim_box_info[5],
&spmd_handle->system->box );
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
x[0] = qm_pos[3 * i];
x[1] = qm_pos[3 * i + 1];
x[2] = qm_pos[3 * i + 2];
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';
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 );
spmd_handle->system->atoms[i].q = 0.0;
spmd_handle->system->atoms[i].qmmm_mask = TRUE;
}
for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
{
x[0] = mm_pos[3 * (i - spmd_handle->system->N_qm)];
x[1] = mm_pos[3 * (i - spmd_handle->system->N_qm) + 1];
x[2] = mm_pos[3 * (i - spmd_handle->system->N_qm) + 2];
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';
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 );
spmd_handle->system->atoms[i].q = mm_q[i - spmd_handle->system->N_qm];
spmd_handle->system->atoms[i].qmmm_mask = FALSE;
}
Read_Input_Files( NULL, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
if ( spmd_handle->system->N > spmd_handle->system->N_max )
{
/* deallocate everything which needs more space
* (i.e., structures whose space is a function of the number of atoms),
* except for data structures allocated while parsing input files */
Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->output_enabled, TRUE );
spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
spmd_handle->realloc = TRUE;
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for atom positions in QMMM mode
*
* handle: pointer to wrapper struct with top-level data structures
* qm_pos: coordinates of QM atom positions (consecutively arranged), in Angstroms (allocated by caller)
* mm_pos: coordinates of MM atom positions (consecutively arranged), in Angstroms (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_positions_qmmm( const void * const handle, double * const qm_pos,
double * const mm_pos )
{
int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
if ( qm_pos != NULL )
{
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
qm_pos[3 * i] = spmd_handle->system->atoms[i].x[0];
qm_pos[3 * i + 1] = spmd_handle->system->atoms[i].x[1];
qm_pos[3 * i + 2] = spmd_handle->system->atoms[i].x[2];
}
}
if ( mm_pos != NULL )
{
for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
{
mm_pos[3 * (i - spmd_handle->system->N_qm)] = spmd_handle->system->atoms[i].x[0];
mm_pos[3 * (i - spmd_handle->system->N_qm) + 1] = spmd_handle->system->atoms[i].x[1]; mm_pos[3 * (i - spmd_handle->system->N_qm) + 2] = spmd_handle->system->atoms[i].x[2];
}
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for atom velocities in QMMM mode
*
* handle: pointer to wrapper struct with top-level data structures
* qm_vel: coordinates of QM atom velocities (consecutively arranged), in Angstroms / ps (allocated by caller)
* mm_vel: coordinates of MM atom velocities (consecutively arranged), in Angstroms / ps (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_velocities_qmmm( const void * const handle, double * const qm_vel,
double * const mm_vel )
{
int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
if ( qm_vel != NULL )
{
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
qm_vel[3 * i] = spmd_handle->system->atoms[i].v[0];
qm_vel[3 * i + 1] = spmd_handle->system->atoms[i].v[1];
qm_vel[3 * i + 2] = spmd_handle->system->atoms[i].v[2];
}
}
if ( mm_vel != NULL )
{
for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
{
mm_vel[3 * (i - spmd_handle->system->N_qm)] = spmd_handle->system->atoms[i].v[0];
mm_vel[3 * (i - spmd_handle->system->N_qm) + 1] = spmd_handle->system->atoms[i].v[1];
mm_vel[3 * (i - spmd_handle->system->N_qm) + 2] = spmd_handle->system->atoms[i].v[2];
}
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for atom forces in QMMM mode
*
* handle: pointer to wrapper struct with top-level data structures
* qm_f: coordinates of QM atom forces (consecutively arranged), in Angstroms * Daltons / ps^2 (allocated by caller)
* mm_f: coordinates of MM atom forces (consecutively arranged), in Angstroms * Daltons / ps^2 (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_forces_qmmm( const void * const handle, double * const qm_f,
double * const mm_f )
{
int i, ret; spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
if ( qm_f != NULL )
{
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
qm_f[3 * i] = spmd_handle->system->atoms[i].f[0];
qm_f[3 * i + 1] = spmd_handle->system->atoms[i].f[1];
qm_f[3 * i + 2] = spmd_handle->system->atoms[i].f[2];
}
}
if ( mm_f != NULL )
{
for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
{
mm_f[3 * (i - spmd_handle->system->N_qm)] = spmd_handle->system->atoms[i].f[0];
mm_f[3 * (i - spmd_handle->system->N_qm) + 1] = spmd_handle->system->atoms[i].f[1];
mm_f[3 * (i - spmd_handle->system->N_qm) + 2] = spmd_handle->system->atoms[i].f[2];
}
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
/* Getter for atom charges in QMMM mode
*
* handle: pointer to wrapper struct with top-level data structures
* qm_q: QM atom charges, in Coulombs (allocated by caller)
* mm_q: MM atom charges, in Coulombs (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_charges_qmmm( const void * const handle, double * const qm_q,
double * const mm_q )
{
int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
if ( qm_q != NULL )
{
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
qm_q[i] = spmd_handle->system->atoms[i].q;
}
}
if ( mm_q != NULL )
{
for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
{
mm_q[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].q;
}
}
ret = SPUREMD_SUCCESS;
}
return ret;
}
#endif
#if defined(QMMM_FORTRAN_AMBER)
/* Allocate top-level data structures and parse input files
* for the first simulation
*
* 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_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_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)
* - angles per dimension (3 entries)
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*/
void setup_qmmm_( void ** handle, const int * const qm_num_atoms,
const int * const qm_types, const double * const qm_pos,
const int * const mm_num_atoms, const int * const mm_types,
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];
rvec x;
spuremd_handle *spmd_handle;
/* top-level allocation */
spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
"setup::spmd_handle" );
/* second-level allocations */
spmd_handle->system = smalloc( sizeof(reax_system),
"Setup::spmd_handle->system" );
spmd_handle->system->prealloc_allocated = FALSE;
spmd_handle->system->ffield_params_allocated = FALSE;
spmd_handle->system->g.allocated = FALSE;
spmd_handle->control = smalloc( sizeof(control_params),
"Setup::spmd_handle->control" );
spmd_handle->data = smalloc( sizeof(simulation_data),
"Setup::spmd_handle->data" );
spmd_handle->workspace = smalloc( sizeof(static_storage),
"Setup::spmd_handle->workspace" );
spmd_handle->workspace->H.allocated = FALSE;
spmd_handle->workspace->H_full.allocated = FALSE;
spmd_handle->workspace->H_sp.allocated = FALSE;
spmd_handle->workspace->H_p.allocated = FALSE;
spmd_handle->workspace->H_spar_patt.allocated = FALSE;
spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
spmd_handle->workspace->H_app_inv.allocated = FALSE;
spmd_handle->workspace->L.allocated = FALSE;
spmd_handle->workspace->U.allocated = FALSE;
spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
"Setup::spmd_handle->lists" );
for ( i = 0; i < LIST_N; ++i )
{ spmd_handle->lists[i] = smalloc( sizeof(reax_list),
"Setup::spmd_handle->lists[i]" );
spmd_handle->lists[i]->allocated = FALSE;
}
spmd_handle->out_control = smalloc( sizeof(output_controls),
"Setup::spmd_handle->out_control" );
spmd_handle->output_enabled = FALSE;
spmd_handle->realloc = TRUE;
spmd_handle->callback = NULL;
spmd_handle->data->sim_id = 0;
spmd_handle->system->N_qm = *qm_num_atoms;
spmd_handle->system->N_mm = *mm_num_atoms;
spmd_handle->system->N = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
PreAllocate_Space( spmd_handle->system, spmd_handle->control,
spmd_handle->workspace, spmd_handle->system->N );
Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
sim_box_info[3], sim_box_info[4], sim_box_info[5],
&spmd_handle->system->box );
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
x[0] = qm_pos[3 * i];
x[1] = qm_pos[3 * i + 1];
x[2] = qm_pos[3 * i + 2];
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';
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 );
spmd_handle->system->atoms[i].q = 0.0;
spmd_handle->system->atoms[i].q_init = 0.0;
spmd_handle->system->atoms[i].qmmm_mask = TRUE;
}
for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
{
x[0] = mm_pos_q[4 * (i - spmd_handle->system->N_qm)];
x[1] = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 1];
x[2] = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 2];
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';
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 );
spmd_handle->system->atoms[i].q = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 3];
spmd_handle->system->atoms[i].q_init = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 3];
spmd_handle->system->atoms[i].qmmm_mask = FALSE;
}
Read_Input_Files( NULL, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
*handle = (void *) spmd_handle;
}
/* Reset for the next simulation by parsing input files and triggering
* reallocation if more space is needed
*
* 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_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_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)
* - angles per dimension (3 entries)
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*/
void reset_qmmm_( const void * const handle,
const int * const qm_num_atoms, const int * const qm_types,
const double * const qm_pos, const int * const mm_num_atoms,
const int * const mm_types, 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;
rvec x;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
/* close files used in previous simulation */
if ( spmd_handle->output_enabled == TRUE )
{
Finalize_Out_Controls( spmd_handle->system, spmd_handle->control,
spmd_handle->workspace, spmd_handle->out_control );
}
spmd_handle->realloc = FALSE;
spmd_handle->data->sim_id++;
spmd_handle->system->N_qm = *qm_num_atoms;
spmd_handle->system->N_mm = *mm_num_atoms;
spmd_handle->system->N = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
PreAllocate_Space( spmd_handle->system, spmd_handle->control,
spmd_handle->workspace, spmd_handle->system->N );
Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
sim_box_info[3], sim_box_info[4], sim_box_info[5],
&spmd_handle->system->box );
for ( i = 0; i < spmd_handle->system->N_qm; ++i )
{
x[0] = qm_pos[3 * i];
x[1] = qm_pos[3 * i + 1]; x[2] = qm_pos[3 * i + 2];
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';
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 );
spmd_handle->system->atoms[i].q = 0.0;
spmd_handle->system->atoms[i].qmmm_mask = TRUE;
}
for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
{
x[0] = mm_pos_q[4 * (i - spmd_handle->system->N_qm)];
x[1] = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 1];
x[2] = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 2];
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';
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 );
spmd_handle->system->atoms[i].q = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 3];
spmd_handle->system->atoms[i].qmmm_mask = FALSE;
}
Read_Input_Files( NULL, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
if ( spmd_handle->system->N > spmd_handle->system->N_max )
{
/* deallocate everything which needs more space
* (i.e., structures whose space is a function of the number of atoms),
* except for data structures allocated while parsing input files */
Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->output_enabled, TRUE );
spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
spmd_handle->realloc = TRUE;
}
ret = SPUREMD_SUCCESS;
}
if ( ret != SPUREMD_SUCCESS )
{
/* TODO: pass errors via another mechanism */
;
}
}
/* Run the simulation according to the prescribed parameters
*
* handle: pointer to wrapper struct with top-level data structures
*/
void simulate_( const void * const handle )
{
int ret;
ret = simulate( handle );
if ( ret != SPUREMD_SUCCESS )
{
/* TODO: pass errors via another mechanism */
;
}
}
/* Deallocate all data structures post-simulation
*
* handle: pointer to wrapper struct with top-level data structures
*/
void cleanup_( const void * const handle )
{
int ret;
ret = cleanup( handle );
if ( ret != SPUREMD_SUCCESS )
{
/* TODO: pass errors via another mechanism */
;
}
}
/* Setter for writing output to files
*
* handle: pointer to wrapper struct with top-level data structures
* enabled: TRUE enables writing output to files, FALSE otherwise
*/
void set_output_enabled_( const void * const handle, const int enabled )
{
int ret;
ret = set_output_enabled( handle, enabled );
if ( ret != SPUREMD_SUCCESS )
{
/* TODO: pass errors via another mechanism */
;
}
}
/* Setter for simulation parameter values as defined in the input control file
*
* handle: pointer to wrapper struct with top-level data structures
* control_keyword: keyword from the control file to set the value for
* control_value: value to set
*/
void set_control_parameter_( const void * const handle, const char * const keyword,
const char ** const values )
{
int ret;
ret = set_control_parameter( handle, keyword, values );
if ( ret != SPUREMD_SUCCESS ) {
/* TODO: pass errors via another mechanism */
;
}
}
/* Getter for atom forces in QMMM mode
*
* handle: pointer to wrapper struct with top-level data structures
* qm_f: coordinates of QM atom forces (consecutively arranged), in Angstroms * Daltons / ps^2 (allocated by caller)
* mm_f: coordinates of MM atom forces (consecutively arranged), in Angstroms * Daltons / ps^2 (allocated by caller)
*/
void get_atom_forces_qmmm_( const void * const handle, double * const qm_f,
double * const mm_f )
{
int ret;
ret = get_atom_forces_qmmm( handle, qm_f, mm_f );
if ( ret != SPUREMD_SUCCESS )
{
/* TODO: pass errors via another mechanism */
;
}
}
/* Getter for atom charges in QMMM mode
*
* handle: pointer to wrapper struct with top-level data structures
* qm_q: QM atom charges, in Coulombs (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
void get_atom_charges_qmmm_( const void * const handle, double * const qm_q )
{
int ret;
ret = get_atom_charges_qmmm( handle, qm_q, NULL );
if ( ret != SPUREMD_SUCCESS )
{
/* TODO: pass errors via another mechanism */
;
}
}
/* Getter for system energies
*
* handle: pointer to wrapper struct with top-level data structures
* e_tot: system total energy, in kcal / mol (reference from caller)
*/
void get_system_info_( const void * const handle, double * const e_tot )
{
int ret;
ret = get_system_info( handle, NULL, NULL, e_tot, NULL, NULL, NULL );
if ( ret != SPUREMD_SUCCESS )
{
/* TODO: pass errors via another mechanism */
;
}
}
#endif