forces.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"

#include "index_utils.h"
#ifdef HAVE_CUDA
#include "cuda_forces.h"

#include "cuda_linear_solvers.h"
#include "cuda_neighbors.h"
//#include "cuda_bond_orders.h"
#include "validation.h"
#endif

#if defined(PURE_REAX)
#include "forces.h"
#include "bond_orders.h"
#include "bonds.h"
#include "basic_comm.h"
#include "hydrogen_bonds.h"
#include "io_tools.h"
#include "list.h"
#include "lookup.h"
#include "multi_body.h"
#include "nonbonded.h"
#include "qEq.h"
#include "tool_box.h"
#include "torsion_angles.h"
#include "valence_angles.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reax_forces.h"
#include "reax_bond_orders.h"
#include "reax_bonds.h"
#include "reax_basic_comm.h"
#include "reax_hydrogen_bonds.h"
#include "reax_io_tools.h"
#include "reax_list.h"
#include "reax_lookup.h"
#include "reax_multi_body.h"
#include "reax_nonbonded.h"
#include "reax_tool_box.h"
#include "reax_torsion_angles.h"
#include "reax_valence_angles.h"
#include "reax_vector.h"
#endif


#ifdef HAVE_CUDA
void Cuda_Total_Forces (reax_system *, control_params *, simulation_data *, storage *);
void Cuda_Total_Forces_PURE (reax_system *, storage *);
#endif


interaction_function Interaction_Functions[NUM_INTRS];


void Dummy_Interaction( reax_system *system, control_params *control,
        simulation_data *data, storage *workspace, reax_list **lists,
        output_controls *out_control )
{
}


void Init_Force_Functions( control_params *control )
{
    Interaction_Functions[0] = BO;
    Interaction_Functions[1] = Bonds; //Dummy_Interaction;
    Interaction_Functions[2] = Atom_Energy; //Dummy_Interaction;
    Interaction_Functions[3] = Valence_Angles; //Dummy_Interaction;
    Interaction_Functions[4] = Torsion_Angles; //Dummy_Interaction;
    if ( control->hbond_cut > 0 )
        Interaction_Functions[5] = Hydrogen_Bonds;
    else Interaction_Functions[5] = Dummy_Interaction;
    Interaction_Functions[6] = Dummy_Interaction; //empty
    Interaction_Functions[7] = Dummy_Interaction; //empty
    Interaction_Functions[8] = Dummy_Interaction; //empty
    Interaction_Functions[9] = Dummy_Interaction; //empty
}


void Compute_Bonded_Forces( reax_system *system, control_params *control,
                            simulation_data *data, storage *workspace,
                            reax_list **lists, output_controls *out_control )
{
    int i;

    /* Mark beginning of a new timestep in bonded energy files */
#if defined(TEST_ENERGY)
    Debug_Marker_Bonded( out_control, data->step );
#endif

    /* Implement all force calls as function pointers */
//  for( i = 0; i < NUM_INTRS; i++ ) {
//#if defined(DEBUG)
//    fprintf( stderr, "p%d: starting f%d\n", system->my_rank, i );
//    MPI_Barrier( MPI_COMM_WORLD );
//#endif
//    (Interaction_Functions[i])( system, control, data, workspace,
//              lists, out_control );
//#if defined(DEBUG)
//    fprintf( stderr, "p%d: f%d done\n", system->my_rank, i );
//    MPI_Barrier( MPI_COMM_WORLD );
//#endif
//  }

    (Interaction_Functions[0])( system, control, data, workspace, lists, out_control );
    (Interaction_Functions[1])( system, control, data, workspace, lists, out_control );
    (Interaction_Functions[2])( system, control, data, workspace, lists, out_control );
    (Interaction_Functions[3])( system, control, data, workspace, lists, out_control );
    (Interaction_Functions[4])( system, control, data, workspace, lists, out_control );
    (Interaction_Functions[5])( system, control, data, workspace, lists, out_control );
}


void Compute_NonBonded_Forces( reax_system *system, control_params *control,
                               simulation_data *data, storage *workspace,
                               reax_list **lists, output_controls *out_control,
                               mpi_datatypes *mpi_data )
{
    /* Mark beginning of a new timestep in nonbonded energy files */
#if defined(TEST_ENERGY)
    Debug_Marker_Nonbonded( out_control, data->step );
#endif

    /* van der Waals and Coulomb interactions */
    if ( control->tabulate == 0 )
        vdW_Coulomb_Energy( system, control, data, workspace,
                            lists, out_control );
    else
        Tabulated_vdW_Coulomb_Energy( system, control, data, workspace,
                                      lists, out_control );

#if defined(DEBUG)
    fprintf( stderr, "p%d: nonbonded forces done\n", system->my_rank );
    MPI_Barrier( MPI_COMM_WORLD );
#endif
}



/* this version of Compute_Total_Force computes forces from
   coefficients accumulated by all interaction functions.
   Saves enormous time & space! */
void Compute_Total_Force( reax_system *system, control_params *control,
                          simulation_data *data, storage *workspace,
                          reax_list **lists, mpi_datatypes *mpi_data )
{
    int i, pj;
    reax_list *bonds = (*lists) + BONDS;

    for ( i = 0; i < system->N; ++i )
        for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
            if ( i < bonds->select.bond_list[pj].nbr )
            {
                if ( control->virial == 0 )
                    Add_dBond_to_Forces( i, pj, workspace, lists );
                else
                    Add_dBond_to_Forces_NPT( i, pj, data, workspace, lists );
            }

    //Print_Total_Force( system, data, workspace );
#if defined(PURE_REAX)
    /* now all forces are computed to their partially-final values
       based on the neighbors information each processor has had.
       final values of force on each atom needs to be computed by adding up
       all partially-final pieces */
    Coll( system, mpi_data, workspace->f, mpi_data->mpi_rvec,
          sizeof(rvec) / sizeof(void), rvec_unpacker );
    for ( i = 0; i < system->n; ++i )
        rvec_Copy( system->my_atoms[i].f, workspace->f[i] );

#if defined(TEST_FORCES)
    Coll( system, mpi_data, workspace->f_ele, mpi_data->mpi_rvec, rvec_unpacker);
    Coll( system, mpi_data, workspace->f_vdw, mpi_data->mpi_rvec, rvec_unpacker);
    Coll( system, mpi_data, workspace->f_be, mpi_data->mpi_rvec, rvec_unpacker );
    Coll( system, mpi_data, workspace->f_lp, mpi_data->mpi_rvec, rvec_unpacker );
    Coll( system, mpi_data, workspace->f_ov, mpi_data->mpi_rvec, rvec_unpacker );
    Coll( system, mpi_data, workspace->f_un, mpi_data->mpi_rvec, rvec_unpacker );
    Coll( system, mpi_data, workspace->f_ang, mpi_data->mpi_rvec, rvec_unpacker);
    Coll( system, mpi_data, workspace->f_coa, mpi_data->mpi_rvec, rvec_unpacker);