init_md.c

/*----------------------------------------------------------------------
  PuReMD - Purdue ReaxFF Molecular Dynamics Program

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

#ifdef HAVE_CUDA
#include "dev_alloc.h"
#include "dev_list.h"
#include "cuda_copy.h"
#include "validation.h"
#endif

#if defined(PURE_REAX)
#include "init_md.h"
#include "allocate.h"
#include "box.h"
#include "comm_tools.h"
#include "forces.h"
#include "grid.h"
#include "integrate.h"
#include "io_tools.h"
#include "list.h"
#include "lookup.h"
#include "neighbors.h"
#include "random.h"
#include "reset_tools.h"
#include "system_props.h"
#include "tool_box.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reax_init_md.h"
#include "reax_allocate.h"
#include "reax_forces.h"
#include "reax_io_tools.h"
#include "reax_list.h"
#include "reax_lookup.h"
#include "reax_reset_tools.h"
#include "reax_system_props.h"
#include "reax_tool_box.h"
#include "reax_vector.h"
#endif


#if defined(PURE_REAX)
/************************ initialize system ************************/
int Reposition_Atoms( reax_system *system, control_params *control,
                      simulation_data *data, mpi_datatypes *mpi_data,
                      char *msg )
{
    int   i;
    rvec  dx;

    /* reposition atoms */
    if ( control->reposition_atoms == 0 )  //fit atoms to periodic box
    {
        rvec_MakeZero( dx );
    }
    else if ( control->reposition_atoms == 1 )  //put center of mass to center
    {
        rvec_Scale( dx, 0.5, system->big_box.box_norms );
        rvec_ScaledAdd( dx, -1., data->xcm );
    }
    else if ( control->reposition_atoms == 2 )  //put center of mass to origin
    {
        rvec_Scale( dx, -1., data->xcm );
    }
    else
    {
        strcpy( msg, "reposition_atoms: invalid option" );
        return FAILURE;
    }

    for ( i = 0; i < system->n; ++i )
        // Inc_on_T3_Gen( system->my_atoms[i].x, dx, &(system->big_box) );
        rvec_Add( system->my_atoms[i].x, dx );

    return SUCCESS;
}



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


    if ( T <= 0.1 )
    {
        for ( i = 0; i < system->n; i++ )
            rvec_MakeZero( system->my_atoms[i].v );
    }
    else
    {
        Randomize();

        for ( i = 0; i < system->n; i++ )
        {
            rvec_Random( system->my_atoms[i].v );

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

            rvec_Scale( system->my_atoms[i].v, 1. / scale, system->my_atoms[i].v );

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

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


int Init_System( reax_system *system, control_params *control,
                 simulation_data *data, storage *workspace,
                 mpi_datatypes *mpi_data, char *msg )
{
    int i;
    reax_atom *atom;
    int nrecv[MAX_NBRS];

    Setup_New_Grid( system, control, MPI_COMM_WORLD );
#if defined(DEBUG_FOCUS)
    fprintf( stderr, "p%d GRID:\n", system->my_rank );
    Print_Grid( &(system->my_grid), stderr );
#endif
    Bin_My_Atoms( system, &(workspace->realloc) );
    Reorder_My_Atoms( system, workspace );

    /* estimate N and total capacity */
    for ( i = 0; i < MAX_NBRS; ++i ) nrecv[i] = 0;
    MPI_Barrier( MPI_COMM_WORLD );
    system->N = SendRecv( system, mpi_data, mpi_data->boundary_atom_type, nrecv,
                          Estimate_Boundary_Atoms, Unpack_Estimate_Message, 1 );
    system->total_cap = MAX( (int)(system->N * SAFE_ZONE), MIN_CAP );
    Bin_Boundary_Atoms( system );

    /* estimate numH and Hcap */
    system->numH = 0;
    if ( control->hbond_cut > 0 )
        for ( i = 0; i < system->n; ++i )
        {
            atom = &(system->my_atoms[i]);
            if ( system->reax_param.sbp[ atom->type ].p_hbond == 1 )
                atom->Hindex = system->numH++;
            else atom->Hindex = -1;
        }
    system->Hcap = MAX( system->numH * SAFER_ZONE, MIN_CAP );

    //Allocate_System( system, system->local_cap, system->total_cap, msg );
#if defined(DEBUG_FOCUS)
    fprintf( stderr, "p%d: n=%d local_cap=%d\n",
             system->my_rank, system->n, system->local_cap );
    fprintf( stderr, "p%d: N=%d total_cap=%d\n",
             system->my_rank, system->N, system->total_cap );
    fprintf( stderr, "p%d: numH=%d H_cap=%d\n",
             system->my_rank, system->numH, system->Hcap );
#endif

    Compute_Total_Mass( system, data, mpi_data->comm_mesh3D );
    Compute_Center_of_Mass( system, data, mpi_data, mpi_data->comm_mesh3D );
    // if( Reposition_Atoms( system, control, data, mpi_data, msg ) == FAILURE )
    //   return FAILURE;

    /* 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 );
    Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );

    return SUCCESS;
}


#ifdef HAVE_CUDA