diff --git a/bin/main.cpp b/bin/main.cpp index f5064bee..f0bf72c6 100644 --- a/bin/main.cpp +++ b/bin/main.cpp @@ -60,6 +60,7 @@ int main( int argc, char *argv[] ) Kokkos::ScopeGuard scope_guard( argc, argv ); + { CabanaMD cabanamd; cabanamd.init( argc, argv ); @@ -68,8 +69,7 @@ int main( int argc, char *argv[] ) // cabanamd.check_correctness(); cabanamd.print_performance(); - - cabanamd.shutdown(); + } // cabanamd's destructor called MPI_Finalize(); } diff --git a/src/cabanamd.cpp b/src/cabanamd.cpp index 529e0379..046d4923 100644 --- a/src/cabanamd.cpp +++ b/src/cabanamd.cpp @@ -55,6 +55,10 @@ #define MAXPATHLEN 1024 CabanaMD::CabanaMD() + : integrator( nullptr ) + , force( nullptr ) + , comm( nullptr ) + , binning( nullptr ) { // Create the System class: atom properties (AoSoA) and simulation box system = new System(); @@ -102,10 +106,13 @@ void CabanaMD::init( int argc, char *argv[] ) // Create Binning class: linked cell bin sort binning = new Binning( system ); + // Create Communication class: MPI + comm = new Comm( system, neigh_cutoff ); + // Create Force class: potential options in force_types/ folder if ( false ) { - } + } // NOTE: modules_force cannot change neigh_cutoff! #define FORCE_MODULES_INSTANTIATION #include #undef FORCE_MODULES_INSTANTIATION @@ -119,9 +126,6 @@ void CabanaMD::init( int argc, char *argv[] ) input->input_data.words[input->force_coeff_lines( line )] ); } - // Create Communication class: MPI - comm = new Comm( system, neigh_cutoff ); - force->comm_newton = input->comm_newton; // system->print_particles(); @@ -564,4 +568,17 @@ void CabanaMD::check_correctness( int step ) void CabanaMD::print_performance() {} -void CabanaMD::shutdown() { system->destroy(); } +CabanaMD::~CabanaMD() +{ + delete system; + delete input; + + if ( integrator != nullptr ) + delete integrator; + if ( force != nullptr ) + delete force; + if ( comm != nullptr ) + delete comm; + if ( binning != nullptr ) + delete binning; +} diff --git a/src/cabanamd.h b/src/cabanamd.h index aac7abd4..3f65a05d 100644 --- a/src/cabanamd.h +++ b/src/cabanamd.h @@ -61,13 +61,14 @@ class CabanaMD { public: System *system; + Input *input; Integrator *integrator; Force *force; Comm *comm; - Input *input; Binning *binning; CabanaMD(); + ~CabanaMD(); void init( int argc, char *argv[] ); @@ -77,6 +78,4 @@ class CabanaMD void check_correctness( int ); void print_performance(); - - void shutdown(); }; diff --git a/src/force.h b/src/force.h index 28addb5c..04019c3b 100644 --- a/src/force.h +++ b/src/force.h @@ -57,6 +57,7 @@ class Force public: bool half_neigh, comm_newton; Force( System *system, bool half_neigh_ ); + virtual ~Force() {}; virtual void init_coeff( T_X_FLOAT neigh_cut, char **args ); virtual void create_neigh_list( System *system ); diff --git a/src/input.cpp b/src/input.cpp index c73c74ca..ca462460 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -187,7 +187,7 @@ Input::Input( System *p ) box[5] = 40; units_style = UNITS_LJ; - lattice_style = LATTICE_FCC; + lattice_style = Lattice::FCC; lattice_constant = 0.8442; temperature_target = 1.4; @@ -351,11 +351,12 @@ void Input::read_lammps_file( const char *filename ) printf( "\n" ); } for ( int l = 0; l < input_data.nlines; l++ ) - check_lammps_command( l ); + run_lammps_command( l ); } -void Input::check_lammps_command( int line ) +void Input::run_lammps_command( int line ) { + System &s = *system; bool known = false; if ( input_data.words[line][0][0] == 0 ) @@ -368,7 +369,7 @@ void Input::check_lammps_command( int line ) } if ( strcmp( input_data.words[line][0], "variable" ) == 0 ) { - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: 'variable' keyword is not supported in " "CabanaMD\n" ); } @@ -378,34 +379,34 @@ void Input::check_lammps_command( int line ) { known = true; units_style = UNITS_METAL; - system->boltz = 8.617343e-5; + s.boltz = 8.617343e-5; // hplanck = 95.306976368; - system->mvv2e = 1.0364269e-4; - system->dt = 0.001; + s.mvv2e = 1.0364269e-4; + s.dt = 0.001; } else if ( strcmp( input_data.words[line][1], "real" ) == 0 ) { known = true; units_style = UNITS_REAL; - system->boltz = 0.0019872067; + s.boltz = 0.0019872067; // hplanck = 95.306976368; - system->mvv2e = 48.88821291 * 48.88821291; + s.mvv2e = 48.88821291 * 48.88821291; if ( !timestepflag ) - system->dt = 1.0; + s.dt = 1.0; } else if ( strcmp( input_data.words[line][1], "lj" ) == 0 ) { known = true; units_style = UNITS_LJ; - system->boltz = 1.0; + s.boltz = 1.0; // hplanck = 0.18292026; - system->mvv2e = 1.0; + s.mvv2e = 1.0; if ( !timestepflag ) - system->dt = 0.005; + s.dt = 0.005; } else { - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: 'units' command only supports " "'metal', 'real', and 'lj' in CabanaMD\n" ); } @@ -419,40 +420,47 @@ void Input::check_lammps_command( int line ) else if ( strcmp( input_data.words[line][1], "charge" ) == 0 ) { known = true; - system->atom_style = "charge"; + s.atom_style = "charge"; } else { - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: 'atom_style' command only supports " "'atomic' and 'charge' in CabanaMD\n" ); } } if ( strcmp( input_data.words[line][0], "lattice" ) == 0 ) { + float atoms_per_unit = 1.0; if ( strcmp( input_data.words[line][1], "sc" ) == 0 ) { known = true; - lattice_style = LATTICE_SC; + lattice_style = Lattice::SC; + lattice_constant = atof( input_data.words[line][2] ); + } + else if ( strcmp( input_data.words[line][1], "bcc" ) == 0 ) + { + known = true; + lattice_style = Lattice::BCC; + atoms_per_unit = 2.0; lattice_constant = atof( input_data.words[line][2] ); } else if ( strcmp( input_data.words[line][1], "fcc" ) == 0 ) { known = true; - lattice_style = LATTICE_FCC; - if ( units_style == UNITS_LJ ) - lattice_constant = - std::pow( ( 4.0 / atof( input_data.words[line][2] ) ), - ( 1.0 / 3.0 ) ); - else - lattice_constant = atof( input_data.words[line][2] ); + lattice_style = Lattice::FCC; + atoms_per_unit = 4.0; + lattice_constant = atof( input_data.words[line][2] ); } else { - if ( system->do_print ) - printf( "LAMMPS-Command: 'lattice' command only supports 'sc' " - "and 'fcc' in CabanaMD\n" ); + if ( s.do_print ) + printf( "LAMMPS-Command: 'lattice' command only supports " + "'sc', 'bcc', or 'fcc' in CabanaMD\n" ); } + if ( units_style == UNITS_LJ ) // words[line][2] was actually rho=N/V + lattice_constant = + std::pow( atoms_per_unit / lattice_constant, 1.0 / 3.0 ); if ( strcmp( input_data.words[line][3], "origin" ) == 0 ) { lattice_offset_x = atof( input_data.words[line][4] ); @@ -465,7 +473,6 @@ void Input::check_lammps_command( int line ) if ( strcmp( input_data.words[line][2], "block" ) == 0 ) { known = true; - int box[6]; box[0] = atoi( input_data.words[line][3] ); box[1] = atoi( input_data.words[line][4] ); box[2] = atoi( input_data.words[line][5] ); @@ -473,7 +480,7 @@ void Input::check_lammps_command( int line ) box[4] = atoi( input_data.words[line][7] ); box[5] = atoi( input_data.words[line][8] ); if ( ( box[0] != 0 ) || ( box[2] != 0 ) || ( box[4] != 0 ) ) - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: region only allows for boxes with " "0,0,0 offset in CabanaMD\n" ); lattice_nx = box[1]; @@ -482,7 +489,7 @@ void Input::check_lammps_command( int line ) } else { - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: 'region' command only supports " "'block' option in CabanaMD\n" ); } @@ -490,8 +497,8 @@ void Input::check_lammps_command( int line ) if ( strcmp( input_data.words[line][0], "create_box" ) == 0 ) { known = true; - system->ntypes = atoi( input_data.words[line][1] ); - system->mass = t_mass( "System::mass", system->ntypes ); + s.ntypes = atoi( input_data.words[line][1] ); + s.mass = t_mass( "System::mass", s.ntypes ); } if ( strcmp( input_data.words[line][0], "create_atoms" ) == 0 ) { @@ -501,7 +508,7 @@ void Input::check_lammps_command( int line ) { known = true; int type = atoi( input_data.words[line][1] ) - 1; - Kokkos::View mass_one( system->mass, type ); + Kokkos::View mass_one( s.mass, type ); T_V_FLOAT mass = atof( input_data.words[line][2] ); Kokkos::deep_copy( mass_one, mass ); } @@ -537,7 +544,7 @@ void Input::check_lammps_command( int line ) Kokkos::resize( force_coeff_lines, 1 ); force_coeff_lines( 0 ) = line; } - if ( system->do_print && !known ) + if ( s.do_print && !known ) printf( "LAMMPS-Command: 'pair_style' command only supports " "'lj/cut' and 'nnp' style in CabanaMD\n" ); } @@ -559,13 +566,13 @@ void Input::check_lammps_command( int line ) known = true; if ( strcmp( input_data.words[line][1], "all" ) != 0 ) { - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: 'velocity' command can only be " "applied to 'all' in CabanaMD\n" ); } if ( strcmp( input_data.words[line][2], "create" ) != 0 ) { - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: 'velocity' command can only be used " "with option 'create' in CabanaMD\n" ); } @@ -595,7 +602,7 @@ void Input::check_lammps_command( int line ) } else { - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: 'fix' command only supports 'nve' " "style in CabanaMD\n" ); } @@ -613,7 +620,7 @@ void Input::check_lammps_command( int line ) if ( strcmp( input_data.words[line][0], "timestep" ) == 0 ) { known = true; - system->dt = atof( input_data.words[line][1] ); + s.dt = atof( input_data.words[line][1] ); timestepflag = true; } if ( strcmp( input_data.words[line][0], "newton" ) == 0 ) @@ -629,7 +636,7 @@ void Input::check_lammps_command( int line ) } else { - if ( system->do_print ) + if ( s.do_print ) printf( "LAMMPS-Command: 'newton' must be followed by 'on' or " "'off'\n" ); } @@ -638,7 +645,7 @@ void Input::check_lammps_command( int line ) { known = true; } - if ( !known && system->do_print ) + if ( !known && s.do_print ) { printf( "ERROR: unknown keyword\n" ); input_data.print_line( line ); @@ -647,133 +654,41 @@ void Input::check_lammps_command( int line ) void Input::create_lattice( Comm *comm ) { - - System s = *system; + System &s = *system; t_mass::HostMirror h_mass = Kokkos::create_mirror_view( s.mass ); Kokkos::deep_copy( h_mass, s.mass ); - // Create Simple Cubic Lattice - if ( lattice_style == LATTICE_SC ) + // Create Simple Cubic Lattice Types + if ( lattice_style == Lattice::SC || lattice_style == Lattice::BCC || + lattice_style == Lattice::FCC ) { - system->domain_x = lattice_constant * lattice_nx; - system->domain_y = lattice_constant * lattice_ny; - system->domain_z = lattice_constant * lattice_nz; - system->domain_hi_x = system->domain_x; - system->domain_hi_y = system->domain_y; - system->domain_hi_z = system->domain_z; + s.domain_x = lattice_constant * lattice_nx; + s.domain_y = lattice_constant * lattice_ny; + s.domain_z = lattice_constant * lattice_nz; + s.domain_hi_x = s.domain_x; + s.domain_hi_y = s.domain_y; + s.domain_hi_z = s.domain_z; comm->create_domain_decomposition(); - s = *system; - - T_INT ix_start = s.sub_domain_lo_x / s.domain_x * lattice_nx - 0.5; - T_INT iy_start = s.sub_domain_lo_y / s.domain_y * lattice_ny - 0.5; - T_INT iz_start = s.sub_domain_lo_z / s.domain_z * lattice_nz - 0.5; - - T_INT ix_end = s.sub_domain_hi_x / s.domain_x * lattice_nx + 0.5; - T_INT iy_end = s.sub_domain_hi_y / s.domain_y * lattice_ny + 0.5; - T_INT iz_end = s.sub_domain_hi_z / s.domain_z * lattice_nz + 0.5; - - T_INT n = 0; - - for ( T_INT iz = iz_start; iz <= iz_end; iz++ ) - { - T_FLOAT ztmp = lattice_constant * ( iz + lattice_offset_z ); - for ( T_INT iy = iy_start; iy <= iy_end; iy++ ) - { - T_FLOAT ytmp = lattice_constant * ( iy + lattice_offset_y ); - for ( T_INT ix = ix_start; ix <= ix_end; ix++ ) - { - T_FLOAT xtmp = lattice_constant * ( ix + lattice_offset_x ); - if ( ( xtmp >= s.sub_domain_lo_x ) && - ( ytmp >= s.sub_domain_lo_y ) && - ( ztmp >= s.sub_domain_lo_z ) && - ( xtmp < s.sub_domain_hi_x ) && - ( ytmp < s.sub_domain_hi_y ) && - ( ztmp < s.sub_domain_hi_z ) ) - { - n++; - } - } - } - } - system->N_local = n; - system->N = n; - system->resize( n ); - s = *system; - auto x = Cabana::slice( s.xvf ); - auto v = Cabana::slice( s.xvf ); - auto id = Cabana::slice( s.xvf ); - auto type = Cabana::slice( s.xvf ); - auto q = Cabana::slice( s.xvf ); - n = 0; - for ( T_INT iz = iz_start; iz <= iz_end; iz++ ) + double basis[4][3] = { + {0.0, 0.0, 0.0}, {0.5, 0.5, 0.0}, {0.5, 0.0, 0.5}, {0.0, 0.5, 0.5}}; + int nbasis = 0; // default = no atoms + switch ( lattice_style ) { - T_FLOAT ztmp = lattice_constant * ( iz + lattice_offset_z ); - for ( T_INT iy = iy_start; iy <= iy_end; iy++ ) - { - T_FLOAT ytmp = lattice_constant * ( iy + lattice_offset_y ); - for ( T_INT ix = ix_start; ix <= ix_end; ix++ ) - { - T_FLOAT xtmp = lattice_constant * ( ix + lattice_offset_x ); - if ( ( xtmp >= s.sub_domain_lo_x ) && - ( ytmp >= s.sub_domain_lo_y ) && - ( ztmp >= s.sub_domain_lo_z ) && - ( xtmp < s.sub_domain_hi_x ) && - ( ytmp < s.sub_domain_hi_y ) && - ( ztmp < s.sub_domain_hi_z ) ) - { - x( n, 0 ) = xtmp; - x( n, 1 ) = ytmp; - x( n, 2 ) = ztmp; - type( n ) = rand() % s.ntypes; - id( n ) = n + 1; - n++; - } - } - } + case Lattice::SC: + nbasis = 1; + break; + case Lattice::BCC: + nbasis = 2; + basis[1][2] = 0.5; // {0,0,0} and {0.5,0.5,0.5} + break; + case Lattice::FCC: + nbasis = 4; + break; } - comm->reduce_int( &system->N, 1 ); - - // Make ids unique over all processes - T_INT N_local_offset = n; - comm->scan_int( &N_local_offset, 1 ); - for ( T_INT i = 0; i < n; i++ ) - id( i ) += N_local_offset - n; - - if ( system->do_print ) - printf( "Atoms: %i %i\n", system->N, system->N_local ); - } - - // Create Face Centered Cubic (FCC) Lattice - if ( lattice_style == LATTICE_FCC ) - { - system->domain_x = lattice_constant * lattice_nx; - system->domain_y = lattice_constant * lattice_ny; - system->domain_z = lattice_constant * lattice_nz; - system->domain_hi_x = system->domain_x; - system->domain_hi_y = system->domain_y; - system->domain_hi_z = system->domain_z; - - comm->create_domain_decomposition(); - s = *system; - - double basis[4][3]; - basis[0][0] = 0.0; - basis[0][1] = 0.0; - basis[0][2] = 0.0; - basis[1][0] = 0.5; - basis[1][1] = 0.5; - basis[1][2] = 0.0; - basis[2][0] = 0.5; - basis[2][1] = 0.0; - basis[2][2] = 0.5; - basis[3][0] = 0.0; - basis[3][1] = 0.5; - basis[3][2] = 0.5; - - for ( int i = 0; i < 4; i++ ) + for ( int i = 0; i < nbasis; i++ ) { basis[i][0] += lattice_offset_x; basis[i][1] += lattice_offset_y; @@ -789,14 +704,14 @@ void Input::create_lattice( Comm *comm ) T_INT iz_end = s.sub_domain_hi_z / s.domain_z * lattice_nz + 0.5; T_INT n = 0; - + // Count local atoms for ( T_INT iz = iz_start; iz <= iz_end; iz++ ) { for ( T_INT iy = iy_start; iy <= iy_end; iy++ ) { for ( T_INT ix = ix_start; ix <= ix_end; ix++ ) { - for ( int k = 0; k < 4; k++ ) + for ( int k = 0; k < nbasis; k++ ) { T_FLOAT xtmp = lattice_constant * ( 1.0 * ix + basis[k][0] ); @@ -817,26 +732,25 @@ void Input::create_lattice( Comm *comm ) } } } + s.N_local = n; + s.N = n; + s.resize( n ); // Allocate space for n local atoms. - system->N_local = n; - system->N = n; - system->resize( n ); - s = *system; auto x = Cabana::slice( s.xvf ); auto v = Cabana::slice( s.xvf ); auto id = Cabana::slice( s.xvf ); auto type = Cabana::slice( s.xvf ); auto q = Cabana::slice( s.xvf ); + // Loop again to fill-in atom properties. n = 0; - for ( T_INT iz = iz_start; iz <= iz_end; iz++ ) { for ( T_INT iy = iy_start; iy <= iy_end; iy++ ) { for ( T_INT ix = ix_start; ix <= ix_end; ix++ ) { - for ( int k = 0; k < 4; k++ ) + for ( int k = 0; k < nbasis; k++ ) { T_FLOAT xtmp = lattice_constant * ( 1.0 * ix + basis[k][0] ); @@ -862,6 +776,7 @@ void Input::create_lattice( Comm *comm ) } } } + comm->reduce_int( &s.N, 1 ); // Make ids unique over all processes T_INT N_local_offset = n; @@ -869,18 +784,16 @@ void Input::create_lattice( Comm *comm ) for ( T_INT i = 0; i < n; i++ ) id( i ) += N_local_offset - n; - comm->reduce_int( &system->N, 1 ); - - if ( system->do_print ) - printf( "Atoms: %i %i\n", system->N, system->N_local ); + if ( s.do_print ) + printf( "Atoms: %i %i\n", s.N, s.N_local ); } + // Initialize velocity using the equivalent of the LAMMPS // velocity geom option, i.e. uniform random kinetic energies. // zero out momentum of the whole system afterwards, to eliminate // drift (bad for energy statistics) - { // Scope s - System s = *system; + { auto x = Cabana::slice( s.xvf ); auto v = Cabana::slice( s.xvf ); auto type = Cabana::slice( s.xvf ); @@ -891,7 +804,7 @@ void Input::create_lattice( Comm *comm ) T_FLOAT total_momentum_y = 0.0; T_FLOAT total_momentum_z = 0.0; - for ( T_INT i = 0; i < system->N_local; i++ ) + for ( T_INT i = 0; i < s.N_local; i++ ) { LAMMPS_RandomVelocityGeom random; double x_i[3] = {x( i, 0 ), x( i, 1 ), x( i, 2 )}; @@ -922,7 +835,7 @@ void Input::create_lattice( Comm *comm ) T_FLOAT system_vy = total_momentum_y / total_mass; T_FLOAT system_vz = total_momentum_z / total_mass; - for ( T_INT i = 0; i < system->N_local; i++ ) + for ( T_INT i = 0; i < s.N_local; i++ ) { v( i, 0 ) -= system_vx; v( i, 1 ) -= system_vy; @@ -934,7 +847,7 @@ void Input::create_lattice( Comm *comm ) T_V_FLOAT T_init_scale = sqrt( temperature_target / T ); - for ( T_INT i = 0; i < system->N_local; i++ ) + for ( T_INT i = 0; i < s.N_local; i++ ) { v( i, 0 ) *= T_init_scale; v( i, 1 ) *= T_init_scale; diff --git a/src/input.h b/src/input.h index 031710d4..f67d0290 100644 --- a/src/input.h +++ b/src/input.h @@ -55,6 +55,14 @@ #include +// Lattice Type +enum class Lattice +{ + SC, + BCC, + FCC +}; + class ItemizedFile { public: @@ -166,7 +174,7 @@ class Input System *system; int units_style; - int lattice_style; + Lattice lattice_style; double lattice_constant, lattice_offset_x, lattice_offset_y, lattice_offset_z; int lattice_nx, lattice_ny, lattice_nz; @@ -215,6 +223,6 @@ class Input void read_file( const char *filename = NULL ); void read_lammps_file( const char *filename ); void read_data_file( const char *filename ); - void check_lammps_command( int line ); + void run_lammps_command( int line ); void create_lattice( Comm *comm ); }; diff --git a/src/integrator_nve.h b/src/integrator_nve.h index aa81cead..58040dd6 100644 --- a/src/integrator_nve.h +++ b/src/integrator_nve.h @@ -63,7 +63,7 @@ class Integrator System *system; Integrator( System *s ); - ~Integrator(); + ~Integrator(){}; T_V_FLOAT timestep_size; void initial_integrate(); diff --git a/src/system.cpp b/src/system.cpp index 09a22588..6b2a69e7 100644 --- a/src/system.cpp +++ b/src/system.cpp @@ -80,7 +80,7 @@ void System::init() mass = t_mass( "System::mass", ntypes ); } -void System::destroy() +System::~System() { N_max = 0; N_local = 0; diff --git a/src/system.h b/src/system.h index 0ed19d68..4572c8f7 100644 --- a/src/system.h +++ b/src/system.h @@ -91,9 +91,8 @@ class System bool print_lammps; System(); - ~System(){}; + ~System(); void init(); - void destroy(); void resize( T_INT new_N ); diff --git a/src/types.h b/src/types.h index bd158790..cd082630 100644 --- a/src/types.h +++ b/src/types.h @@ -62,12 +62,6 @@ enum UNITS_LJ, UNITS_METAL }; -// Lattice Type -enum -{ - LATTICE_SC, - LATTICE_FCC -}; // Integrator Type enum { diff --git a/unit_test/lattice_neighbors.cpp b/unit_test/lattice_neighbors.cpp new file mode 100644 index 00000000..5057f736 --- /dev/null +++ b/unit_test/lattice_neighbors.cpp @@ -0,0 +1,234 @@ +/**************************************************************************** + * Copyright (c) 2018-2019 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +//************************************************************************ +// ExaMiniMD v. 1.0 +// Copyright (2018) National Technology & Engineering Solutions of Sandia, +// LLC (NTESS). +// +// Under the terms of Contract DE-NA-0003525 with NTESS, the U.S. Government +// retains certain rights in this software. +// +// ExaMiniMD is licensed under 3-clause BSD terms of use: Redistribution and +// use in source and binary forms, with or without modification, are +// permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY NTESS "AS IS" AND ANY EXPRESS OR IMPLIED +// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL NTESS OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING +// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +//************************************************************************ + +#include +#include +#include "mpi.h" +#include + +static constexpr T_FLOAT sqrt2 = std::sqrt(2.0); +static constexpr T_FLOAT sqrt3 = std::sqrt(3.0); + +// Count the number of nearest and next-nearest neighbors in a given lattice. +struct Nbrs { + char lattice[4]; + int n1, n2; + T_FLOAT d1, d2; + T_FLOAT atoms_per_unit; // dens = atoms_per_unit / lattice_const^3 +} nbrs[] = { {"sc", 6, 12, 1.0, sqrt2, 1.0}, + {"bcc", 8, 6, sqrt3/2.0, 1.0, 2.0}, + {"fcc", 12, 6, 1.0/sqrt2, 1.0, 4.0}, + {"hcp", 12, 0, 1.0, 1.0, 1.0} + }; + +// E_LJ = 4*eps*( (sigma/r)^12 - (sigma/r)^6 ) + +T_FLOAT compute_pe(CabanaMD &cmd) { + PotE pote( cmd.comm ); + T_FLOAT PE = pote.compute( cmd.system, cmd.force ) / cmd.system->N; + return PE; +} + +void init(CabanaMD &cmd, Lattice lat, + T_FLOAT lattice_constant, T_FLOAT force_cutoff) +{ + // These are needed because modules_force is macro-expanded. + Input *input = cmd.input; + Force * &force = cmd.force; + System *system = cmd.system; + System &s = *system; + + // Parse command line arguments + // Read input file + auto units_style = UNITS_LJ; + + s.boltz = 1.0; + s.mvv2e = 1.0; + s.dt = 0.005; + input->lattice_style = lat; + input->lattice_constant = lattice_constant; + + input->box[1] = 10; + input->box[3] = 10; + input->box[5] = 10; + input->lattice_nx = input->box[1]; + input->lattice_ny = input->box[3]; + input->lattice_nz = input->box[5]; + s.ntypes = 1; + s.mass = t_mass( "System::mass", s.ntypes ); + //input->neighbor_skin = 0.3; + input->neighbor_skin = 0.0; + //input->force_cutoff = force_cutoff; // only used in initialization + + // CabanaMD.init() resumes + T_X_FLOAT neigh_cutoff = force_cutoff + input->neighbor_skin; + cmd.integrator = new Integrator( system ); // NVE integrator + cmd.binning = new Binning( system ); // linked cell bin sort + cmd.comm = new Comm( system, neigh_cutoff ); // MPI Communicator + + if( false ) {} +#define FORCE_MODULES_INSTANTIATION +#include +#undef FORCE_MODULES_INSTANTIATION + else + cmd.comm->error( "Invalid ForceType" ); + + /*for ( std::size_t line = 0; line < input->force_coeff_lines.extent( 0 ); + line++ ) + { + cmd.force->init_coeff( + neigh_cutoff, + input->input_data.words[input->force_coeff_lines( line )] ); + }*/ + char _args[][8] = {"", "", "", + "1.0", // eps + "1.0", // sigma + "3.0", // cut + }; + char *args[6] = {_args[0], _args[1], _args[2], + _args[3], _args[4], _args[5]}; + snprintf(args[5], 8, "%7.4f", force_cutoff); + cmd.force->init_coeff(neigh_cutoff, args); + cmd.force->comm_newton = input->comm_newton; // Newton's 2nd law setting + + // ** create lattice ** + input->create_lattice( cmd.comm ); + + // exchange atoms across all MPI ranks + cmd.comm->exchange(); + + // Sort atoms + cmd.binning->create_binning( neigh_cutoff, neigh_cutoff, neigh_cutoff, + 1, true, false, true ); + + // Add ghost atoms from other MPI ranks (gather) + cmd.comm->exchange_halo(); + + // Compute atom neighbors + cmd.force->create_neigh_list( cmd.system ); + + // Compute initial forces + auto f = Cabana::slice( s.xvf ); + Cabana::deep_copy( f, 0.0 ); + cmd.force->compute( cmd.system ); + + // Scatter ghost atom forces back to original MPI rank + // (update force for pair_style nnp even if full neighbor list) + if ( input->comm_newton or input->force_type == FORCE_NNP ) { + cmd.comm->update_force(); + } +} + +Lattice parse_lattice(const char *style) { + if ( strcmp( style, "fcc" ) == 0 ) + { + return Lattice::FCC; + } + else if ( strcmp( style, "bcc" ) == 0 ) + { + return Lattice::BCC; + } + return Lattice::SC; +} + +T_FLOAT lattice_energy(const char *style, T_FLOAT a, T_FLOAT cut) { + struct Nbrs *nbr = nbrs; + T_FLOAT en; + + for(int i=0; ilattice, 4)) { + break; + } + } + if(nbr >= nbrs+sizeof(nbrs)/sizeof(nbrs[0])) { + return 0.0; + } + if(a*nbr->d1 < cut) { + T_FLOAT ir6 = std::pow(a*nbr->d1, -6.0); + en += nbr->n1*4.0*(ir6*ir6 - ir6); + } + if(a*nbr->d2 < cut) { + T_FLOAT ir6 = std::pow(a*nbr->d2, -6.0); + en += nbr->n2*4.0*(ir6*ir6 - ir6); + } + + return en; +} + +int main( int argc, char *argv[] ) { + MPI_Init( &argc, &argv ); + Kokkos::ScopeGuard scope_guard( argc, argv ); + if(argc != 4) { + printf("Usage: %s \n", argv[0]); + MPI_Finalize(); + return 1; + } + + { + Lattice shape = parse_lattice(argv[1]); + T_FLOAT a = atof(argv[2]); + T_FLOAT cut = atof(argv[3]); + CabanaMD cabanamd; + //cabanamd.init( argc, argv ); + init(cabanamd, shape, a, cut); + + T_FLOAT PE = compute_pe(cabanamd); + printf("Energy = %g\n", PE); + printf("Lattice Energy = %g per atom\n", + lattice_energy(argv[1], a, cut)); + + //cabanamd.run( cabanamd.input->nsteps ); + //cabanamd.check_correctness(); + //cabanamd.print_performance(); + } + + MPI_Finalize(); + + return 0; +} +