spuremd.c

/*----------------------------------------------------------------------
  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 );
    }

    if ( control_file != NULL )
    {
        Read_Control_File( control_file, system, control, out_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
 *
 * qm_num_atoms: num. atoms in the QM region
 * qm_types: element types for QM atoms
 * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
 * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
 * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
 * mm_num_atoms: num. atoms in the MM region
 * mm_types: element types for MM atoms
 * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
 * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
 * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
 * mm_q: charge of MM atom, in Coulombs
 * sim_box: 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_x, const double * const qm_pos_y,
        const double * const qm_pos_z,
        int mm_num_atoms, const int * const mm_types,
        const double * const mm_pos_x, const double * const mm_pos_y,
        const double * const mm_pos_z, const double * const mm_q,
        const double * const sim_box,
        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" );