init_md.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 "init_md.h"

#include "allocate.h"
#include "box.h"
#include "forces.h"
#include "grid.h"
#include "integrate.h"
#include "neighbors.h"
#include "list.h"
#include "lookup.h"
#include "reset_utils.h"
#include "system_props.h"
#include "tool_box.h"
#include "vector.h"


void Generate_Initial_Velocities( reax_system *system, real T )
{
    int i;
    real scale, norm;

    if ( T <= 0.1 )
    {
        for ( i = 0; i < system->N; i++ )
        {
            rvec_MakeZero( system->atoms[i].v );
        }

#if defined(DEBUG)
        fprintf( stderr, "no random velocities...\n" );
#endif
    }
    else
    {
        for ( i = 0; i < system->N; i++ )
        {
            rvec_Random( system->atoms[i].v );

            norm = rvec_Norm_Sqr( system->atoms[i].v );
            scale = SQRT( system->reaxprm.sbp[ system->atoms[i].type ].mass *
                          norm / (3.0 * K_B * T) );

            rvec_Scale( system->atoms[i].v, 1.0 / scale, system->atoms[i].v );

            /*fprintf( stderr, "v = %f %f %f\n",
            system->atoms[i].v[0],system->atoms[i].v[1],system->atoms[i].v[2]);
            fprintf( stderr, "scale = %f\n", scale );
            fprintf( stderr, "v = %f %f %f\n",
            system->atoms[i].v[0],system->atoms[i].v[1],system->atoms[i].v[2]);*/
        }
    }
}


void Init_System( reax_system *system, control_params *control,
        simulation_data *data )
{
    int i;
    rvec dx;

    if ( !control->restart )
    {
        Reset_Atoms( system );
    }

    Compute_Total_Mass( system, data );
    Compute_Center_of_Mass( system, data, stderr );

    /* reposition atoms */
    // just fit the atoms to the periodic box
    if ( control->reposition_atoms == 0 )
    {
        rvec_MakeZero( dx );
    }
    // put the center of mass to the center of the box
    else if ( control->reposition_atoms == 1 )
    {
        rvec_Scale( dx, 0.5, system->box.box_norms );
        rvec_ScaledAdd( dx, -1., data->xcm );
    }
    // put the center of mass to the origin
    else if ( control->reposition_atoms == 2 )
    {
        rvec_Scale( dx, -1., data->xcm );
    }
    else
    {
        fprintf( stderr, "UNKNOWN OPTION: reposition_atoms. Terminating...\n" );
        exit( UNKNOWN_OPTION );
    }

    for ( i = 0; i < system->N; ++i )
    {
        Inc_on_T3( system->atoms[i].x, dx, &(system->box) );
        /*fprintf( stderr, "%6d%2d%8.3f%8.3f%8.3f\n",
          i, system->atoms[i].type,
          system->atoms[i].x[0], system->atoms[i].x[1], system->atoms[i].x[2] );*/
    }

    /* Initialize velocities so that desired init T can be attained */
    if ( !control->restart || (control->restart && control->random_vel) )
    {
        Generate_Initial_Velocities( system, control->T_init );
    }

    Setup_Grid( system );
}


void Init_Simulation_Data( reax_system *system, control_params *control,
        simulation_data *data, output_controls *out_control,
        evolve_function *Evolve )
{

    Reset_Simulation_Data( data );

    if ( !control->restart )
    {
        data->step = 0;
        data->prev_steps = 0;
    }

    switch ( control->ensemble )
    {
    case NVE:
        data->N_f = 3 * system->N;
        *Evolve = Velocity_Verlet_NVE;
        break;


    case NVT:
        data->N_f = 3 * system->N + 1;
        //control->Tau_T = 100 * data->N_f * K_B * control->T_final;
        if ( !control->restart || (control->restart && control->random_vel) )
        {
            data->therm.G_xi = control->Tau_T * (2.0 * data->E_Kin -
                                                 data->N_f * K_B * control->T );
            data->therm.v_xi = data->therm.G_xi * control->dt;
            data->therm.v_xi_old = 0;
            data->therm.xi = 0;
#if defined(DEBUG_FOCUS)
            fprintf( stderr, "init_md: G_xi=%f Tau_T=%f E_kin=%f N_f=%f v_xi=%f\n",
                     data->therm.G_xi, control->Tau_T, data->E_Kin,
                     data->N_f, data->therm.v_xi );
#endif
        }

        *Evolve = Velocity_Verlet_Nose_Hoover_NVT_Klein;
        break;


    case NPT: // Anisotropic NPT
        fprintf( stderr, "THIS OPTION IS NOT YET IMPLEMENTED! TERMINATING...\n" );
        exit( UNKNOWN_OPTION );
        data->N_f = 3 * system->N + 9;
        if ( !control->restart )
        {
            data->therm.G_xi = control->Tau_T * (2.0 * data->E_Kin -
                    data->N_f * K_B * control->T);
            data->therm.v_xi = data->therm.G_xi * control->dt;
            data->iso_bar.eps = 1.0 / 3.0 * LOG( system->box.volume );
            //data->inv_W = 1. / (data->N_f*K_B*control->T*SQR(control->Tau_P));
            //Compute_Pressure( system, data, workspace );
        }
        *Evolve = Velocity_Verlet_Berendsen_Isotropic_NPT;
        break;


    case sNPT: // Semi-Isotropic NPT
        data->N_f = 3 * system->N + 4;
        *Evolve = Velocity_Verlet_Berendsen_SemiIsotropic_NPT;
        break;


    case iNPT: // Isotropic NPT
        data->N_f = 3 * system->N + 2;
        *Evolve = Velocity_Verlet_Berendsen_Isotropic_NPT;
        break;