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
int Cuda_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];