include 2 use cases in doc page

This commit is contained in:
Steve Plimpton
2021-08-18 11:55:47 -06:00
4331 changed files with 85958 additions and 74705 deletions

174
src/.gitignore vendored
View File

@ -52,7 +52,7 @@
/*_dielectric.cpp
/*_dielectric.h
/compute_efield_atom.cpp
/compute_efield_atom.j
/compute_efield_atom.h
/superpose3d.h
@ -70,13 +70,14 @@
/meam*.h
/meam*.cpp
/pair_meamc.cpp
/pair_meamc.h
/pair_meam.cpp
/pair_meam.h
/compute_mliap.cpp
/compute_mliap.h
/mliap_*.cpp
/mliap_*.h
/mliap_model_python_couple.pyx
/pair_mliap.cpp
/pair_mliap.h
@ -250,6 +251,68 @@
/pair_mesont_tpm.cpp
/pair_mesont_tpm.h
/compute_adf.cpp
/compute_adf.h
/compute_contact_atom.cpp
/compute_contact_atom.h
/compute_hexorder_atom.cpp
/compute_hexorder_atom.h
/dump_dcd.cpp
/dump_dcd.h
/fix_controller.cpp
/fix_controller.h
/fix_drag.cpp
/fix_drag.h
/fix_numdiff.cpp
/fix_numdiff.h
/fix_nve_noforce.cpp
/fix_nve_noforce.h
/fix_spring_rg.cpp
/fix_spring_rg.h
/fix_temp_csld.cpp
/fix_temp_csld.h
/fix_temp_csvr.cpp
/fix_temp_csvr.h
/fix_tmd.cpp
/fix_tmd.h
/pair_beck.cpp
/pair_beck.h
/pair_born_coul_dsf.cpp
/pair_born_coul_dsf.h
/pair_born_coul_wolf.cpp
/pair_born_coul_wolf.h
/pair_coul_cut_global.cpp
/pair_coul_cut_global.h
/pair_coul_streitz.cpp
/pair_coul_streitz.h
/pair_gauss.cpp
/pair_gauss.h
/pair_lj96_cut.cpp
/pair_lj96_cut.h
/pair_lj_cubic.cpp
/pair_lj_cubic.h
/pair_lj_cubic_const.h
/pair_lj_cut_coul_debye.cpp
/pair_lj_cut_coul_debye.h
/pair_lj_cut_coul_dsf.cpp
/pair_lj_cut_coul_dsf.h
/pair_lj_cut_coul_wolf.cpp
/pair_lj_cut_coul_wolf.h
/pair_lj_gromacs.cpp
/pair_lj_gromacs.h
/pair_lj_gromacs_coul_gromacs.cpp
/pair_lj_gromacs_coul_gromacs.h
/pair_lj_relres.cpp
/pair_lj_relres.h
/pair_lj_smooth.cpp
/pair_lj_smooth.h
/pair_lj_smooth_linear.cpp
/pair_lj_smooth_linear.h
/pair_mie_cut.cpp
/pair_mie_cut.h
/pair_ufm.cpp
/pair_ufm.h
/angle_cg_cmm.cpp
/angle_cg_cmm.h
/angle_charmm.cpp
@ -314,6 +377,8 @@
/atom_vec_granular.h
/atom_vec_molecular.cpp
/atom_vec_molecular.h
/atom_vec_oxdna.cpp
/atom_vec_oxdna.h
/atom_vec_peri.cpp
/atom_vec_peri.h
/atom_vec_template.cpp
@ -386,6 +451,8 @@
/compute_erotate_rigid.h
/compute_event_displace.cpp
/compute_event_displace.h
/compute_fabric.cpp
/compute_fabric.h
/compute_fep.cpp
/compute_fep.h
/compute_force_tally.cpp
@ -574,12 +641,8 @@
/fix_client_md.h
/fix_cmap.cpp
/fix_cmap.h
/fix_deposit.cpp
/fix_deposit.h
/fix_dpd_energy.cpp
/fix_dpd_energy.h
/fix_efield.cpp
/fix_efield.h
/fix_electron_stopping.cpp
/fix_electron_stopping.h
/fix_electron_stopping_fit.cpp
@ -588,8 +651,6 @@
/fix_eos_cv.h
/fix_eos_table.cpp
/fix_eos_table.h
/fix_evaporate.cpp
/fix_evaporate.h
/fix_ffl.cpp
/fix_ffl.h
/fix_filter_corotate.cpp
@ -718,6 +779,8 @@
/fix_orient_eco.h
/fix_orient_fcc.cpp
/fix_orient_fcc.h
/fix_pair_tracker.cpp
/fix_pair_tracker.h
/fix_peri_neigh.cpp
/fix_peri_neigh.h
/fix_phonon.cpp
@ -728,18 +791,18 @@
/fix_pour.h
/fix_qeq_comb.cpp
/fix_qeq_comb.h
/fix_qeq_reax.cpp
/fix_qeq_fire.cpp
/fix_qeq_fire.h
/fix_qeq_reax.h
/fix_qeq_reaxff.cpp
/fix_qeq_reaxff.h
/fix_qmmm.cpp
/fix_qmmm.h
/fix_reaxc.cpp
/fix_reaxc.h
/fix_reaxc_bonds.cpp
/fix_reaxc_bonds.h
/fix_reaxc_species.cpp
/fix_reaxc_species.h
/fix_reaxff.cpp
/fix_reaxff.h
/fix_reaxff_bonds.cpp
/fix_reaxff_bonds.h
/fix_reaxff_species.cpp
/fix_reaxff_species.h
/fix_rhok.cpp
/fix_rhok.h
/fix_rigid.cpp
@ -791,8 +854,6 @@
/fix_tgnpt_drude.h
/fix_tgnvt_drude.cpp
/fix_tgnvt_drude.h
/fix_thermal_conductivity.cpp
/fix_thermal_conductivity.h
/fix_ti_rs.cpp
/fix_ti_rs.h
/fix_ti_spring.cpp
@ -1129,8 +1190,8 @@
/pair_peri_ves.h
/pair_quip.cpp
/pair_quip.h
/pair_reaxc.cpp
/pair_reaxc.h
/pair_reaxff.cpp
/pair_reaxff.h
/pair_rebo.cpp
/pair_rebo.h
/pair_resquared.cpp
@ -1167,6 +1228,8 @@
/pair_tip4p_long.h
/pair_tip4p_long_soft.cpp
/pair_tip4p_long_soft.h
/pair_tracker.h
/pair_tracker.cpp
/pair_tri_lj.cpp
/pair_tri_lj.h
/pair_yukawa_colloid.cpp
@ -1215,50 +1278,27 @@
/reader_adios.h
/reader_molfile.cpp
/reader_molfile.h
/reaxc_allocate.cpp
/reaxc_allocate.h
/reaxc_basic_comm.cpp
/reaxc_basic_comm.h
/reaxc_bond_orders.cpp
/reaxc_bond_orders.h
/reaxc_bonds.cpp
/reaxc_bonds.h
/reaxc_control.cpp
/reaxc_control.h
/reaxc_defs.h
/reaxc_ffield.cpp
/reaxc_ffield.h
/reaxc_forces.cpp
/reaxc_forces.h
/reaxc_hydrogen_bonds.cpp
/reaxc_hydrogen_bonds.h
/reaxc_init_md.cpp
/reaxc_init_md.h
/reaxc_io_tools.cpp
/reaxc_io_tools.h
/reaxc_list.cpp
/reaxc_list.h
/reaxc_lookup.cpp
/reaxc_lookup.h
/reaxc_multi_body.cpp
/reaxc_multi_body.h
/reaxc_nonbonded.cpp
/reaxc_nonbonded.h
/reaxc_reset_tools.cpp
/reaxc_reset_tools.h
/reaxc_system_props.cpp
/reaxc_system_props.h
/reaxc_tool_box.cpp
/reaxc_tool_box.h
/reaxc_torsion_angles.cpp
/reaxc_torsion_angles.h
/reaxc_traj.cpp
/reaxc_traj.h
/reaxc_types.h
/reaxc_valence_angles.cpp
/reaxc_valence_angles.h
/reaxc_vector.cpp
/reaxc_vector.h
/reaxff_allocate.cpp
/reaxff_bond_orders.cpp
/reaxff_bonds.cpp
/reaxff_control.cpp
/reaxff_ffield.cpp
/reaxff_forces.cpp
/reaxff_hydrogen_bonds.cpp
/reaxff_init_md.cpp
/reaxff_list.cpp
/reaxff_lookup.cpp
/reaxff_multi_body.cpp
/reaxff_nonbonded.cpp
/reaxff_reset_tools.cpp
/reaxff_tool_box.cpp
/reaxff_torsion_angles.cpp
/reaxff_valence_angles.cpp
/reaxff_api.h
/reaxff_defs.h
/reaxff_inline.h
/reaxff_omp.h
/reaxff_types.h
/remap.cpp
/remap.h
/remap_wrap.cpp
@ -1429,6 +1469,10 @@
/pair_thole.h
/pair_buck_mdf.cpp
/pair_buck_mdf.h
/pair_dpd.cpp
/pair_dpd.h
/pair_dpd_tstat.cpp
/pair_dpd_tstat.h
/pair_dpd_ext.cpp
/pair_dpd_ext.h
/pair_dpd_ext_tstat.cpp

View File

@ -72,7 +72,7 @@ adios_SYSLIB=${ADIOS2_LIB}
sed -i -e '/^include.*ADIOS.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ../USER-ADIOS/Makefile.lammps
include ../ADIOS/Makefile.lammps
' ../Makefile.package.settings
fi
fi

View File

@ -8,7 +8,7 @@ Configure LAMMPS with CMake
ADIOS2_DIR
to the ADIOS 2.x installation path
b. use the cmake option
-D PKG_USER-ADIOS=yes
-D PKG_ADIOS=yes
The person who created this package is Norbert Podhorszki (Oak Ridge National Laboratory);
If you need help, please submit a ticket at the OLCF ticket user support mentioning his name in the ticket.

File diff suppressed because it is too large Load Diff

View File

@ -31,7 +31,7 @@ action () {
if (test $1 = 1) then
if (test ! -e ../pair_eam.cpp) then
echo "Must install MANYBODY package with USER-ATC"
echo "Must install MANYBODY package with ATC"
exit 1
fi
fi
@ -54,9 +54,9 @@ if (test $1 = 1) then
# sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/atc |' ../Makefile.package
# sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/atc$(LIBSOBJDIR) |' ../Makefile.package
sed -i -e 's|^PKG_LIB =[ \t]*|&-latc |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(user-atc_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(user-atc_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(user-atc_SYSPATH) |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(atc_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(atc_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(atc_SYSPATH) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then

View File

@ -7,7 +7,7 @@ See the doc page for the fix atc command to get started. At the
bottom of the doc page are many links to additional documentation
contained in the doc/USER/atc directory.
There are example scripts for using this package in examples/USER/atc.
There are example scripts for using this package in examples/PACKAGES/atc.
This package uses an external library in lib/atc which must be
compiled before making LAMMPS. See the lib/atc/README file and the

View File

@ -37,7 +37,7 @@ using namespace FixConst;
using std::string;
#ifdef LAMMPS_BIGBIG
#error "The USER-ATC package is not compatible with -DLAMMPS_BIGBIG"
#error "The ATC package is not compatible with -DLAMMPS_BIGBIG"
#endif
// main page of doxygen documentation
@ -700,20 +700,17 @@ int FixATC::size_restart(int /* nlocal */) {
void FixATC::write_restart(FILE * /* fp */) {
char ** args = new char*[2];
args[0] = new char[50];
args[1] = new char[50];
sprintf(args[0],"write_restart");
sprintf(args[1],"ATC.restart");
char *args[2];
args[0] = utils::strdup("write_restart");
args[1] = utils::strdup("ATC.restart");
// Then call all objects I own to write their data
if (comm->me == 0) {
atc_->modify(2,args);
}
delete [] args[0];
delete [] args[1];
delete [] args;
delete[] args[0];
delete[] args[1];
}
/* ----------------------------------------------------------------------
@ -722,20 +719,17 @@ void FixATC::write_restart(FILE * /* fp */) {
void FixATC::restart(char * /* buf */) {
char ** args = new char*[2];
args[0] = new char[50];
args[1] = new char[50];
sprintf(args[0],"read_restart");
sprintf(args[1],"ATC.restart");
char *args[2];
args[0] = utils::strdup("read_restart");
args[1] = utils::strdup("ATC.restart");
// Then call all objects I own to write their data
if (comm->me == 0) {
atc_->modify(2,args);
}
delete [] args[0];
delete [] args[1];
delete [] args;
delete[] args[0];
delete[] args[1];
}
/* ----------------------------------------------------------------------

View File

@ -41,8 +41,8 @@ if (test $1 = 1) then
sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/awpmd/ivutils/include -I../../lib/awpmd/systems/interact |' ../Makefile.package
sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/awpmd |' ../Makefile.package
sed -i -e 's|^PKG_LIB =[ \t]*|&-lawpmd |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(user-awpmd_SYSPATH) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(user-awpmd_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(awpmd_SYSPATH) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(awpmd_SYSLIB) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then

View File

@ -4,7 +4,7 @@ Wave Packet Molecular Dynamics (AWPMD) method.
See the doc page for the pair_style awpmd/cut command to get started.
There are example scripts for using this package in
examples/USER/awpmd.
examples/PACKAGES/awpmd.
This package uses an external library in lib/awpmd which must be
compiled before making LAMMPS. See the lib/awpmd/README file and the

View File

@ -7,7 +7,7 @@ molecular liquids," J. Chem. Phys. 143, 243148 (2015).
doi: 10.1063/1.4937383
The USER-BOCS user package for LAMMPS is part of the BOCS software package:
The BOCS user package for LAMMPS is part of the BOCS software package:
https://github.com/noid-group/BOCS
See the following reference for information about the entire package:
@ -16,4 +16,4 @@ Dunn, NJH; Lebold, KM; DeLyser, MR; Rudzinski, JF; Noid, WG.
"BOCS: Bottom-Up Open-Source Coarse-Graining Software."
J. Phys. Chem. B. 122, 13, 3363-3377 (2018).
Example inputs are in the examples/USER/bocs folder.
Example inputs are in the examples/PACKAGES/bocs folder.

View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------
USER-BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
from The Pennsylvania State University
------------------------------------------------------------------------- */

View File

@ -10,7 +10,7 @@
See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------
USER-BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
from The Pennsylvania State University
------------------------------------------------------------------------- */

View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------
USER-BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
from The Pennsylvania State University
------------------------------------------------------------------------- */
@ -43,7 +43,7 @@ using namespace LAMMPS_NS;
using namespace FixConst;
static const char cite_user_bocs_package[] =
"USER-BOCS package:\n\n"
"BOCS package:\n\n"
"@Article{Dunn2018,\n"
" author = {NJH Dunn, KM Lebold, MR DeLyser, JF Rudzinski, WG Noid},\n"
" title = {BOCS: Bottom-Up Open-Source Coarse-Graining Software},\n"
@ -1100,7 +1100,12 @@ void FixBocs::initial_integrate_respa(int /*vflag*/, int ilevel, int /*iloop*/)
nve_x();
if (pstat_flag) remap();
}
}
/* ---------------------------------------------------------------------- */
void FixBocs::pre_force_respa(int /*vflag*/, int ilevel, int /*iloop*/)
{
// if barostat, redo KSpace coeffs at outermost level,
// since volume has changed

View File

@ -10,7 +10,7 @@
See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------
USER-BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
from The Pennsylvania State University
------------------------------------------------------------------------- */
@ -37,6 +37,7 @@ class FixBocs : public Fix {
virtual void initial_integrate(int);
virtual void final_integrate();
void initial_integrate_respa(int, int, int);
void pre_force_respa(int, int, int);
void final_integrate_respa(int, int);
virtual void pre_exchange();
double compute_scalar();

View File

@ -21,21 +21,21 @@
#include "pair_body_rounded_polygon.h"
#include <cmath>
#include <cstring>
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_body.h"
#include "body_rounded_polygon.h"
#include "comm.h"
#include "force.h"
#include "fix.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "math_extra.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;

View File

@ -23,22 +23,22 @@
#include "pair_body_rounded_polyhedron.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "atom_vec_body.h"
#include "body_rounded_polyhedron.h"
#include "comm.h"
#include "force.h"
#include "fix.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "math_extra.h"
#include "fix.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -1212,10 +1212,7 @@ int PairBodyRoundedPolyhedron::interaction_edge_to_edge(int ibody,
contact_list[num_contacts].unique = 1;
num_contacts++;
}
} else {
}
return interact;
}

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp.
Originally modified from CG-DNA/fix_nve_dotc_langevin.cpp.
Contributing author: Sam Cameron (University of Bristol)
------------------------------------------------------------------------- */

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp.
Originally modified from CG-DNA/fix_nve_dotc_langevin.cpp.
Contributing author: Sam Cameron (University of Bristol)
------------------------------------------------------------------------- */

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp.
Originally modified from CG-DNA/fix_nve_dotc_langevin.cpp.
Contributing author: Sam Cameron (University of Bristol)
------------------------------------------------------------------------- */

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp.
Originally modified from CG-DNA/fix_nve_dotc_langevin.cpp.
Contributing author: Sam Cameron (University of Bristol)
------------------------------------------------------------------------- */

View File

@ -26,8 +26,25 @@ action () {
fi
}
# the CG-DNA package cannot be used without
# the MOLECULE and ASPHERE packages installed.
if (test $1 = 1) then
if (test ! -e ../bond_harmonic.cpp) then
echo "Must install MOLECULE package with CG-DNA"
exit 1
fi
if (test ! -e ../fix_nve_asphere.cpp) then
echo "Must install ASPHERE package with CG-DNA"
exit 1
fi
fi
# list of files with dependcies
action atom_vec_oxdna.cpp
action atom_vec_oxdna.h
action bond_oxdna_fene.cpp bond_fene.h
action bond_oxdna2_fene.cpp bond_fene.h
action bond_oxrna2_fene.cpp bond_fene.h

View File

@ -40,13 +40,13 @@ gradient thermostats for rigid body dynamics", J. Chem. Phys. 142,
144114 (2015).
Example input and data files can be found in
/examples/USER/cgdna/examples/oxDNA/, /oxDNA2/ and /oxRNA2/.
examples/PACKAGES/cgdna/examples/oxDNA/, /oxDNA2/ and /oxRNA2/.
Python setup tools which create single straight or helical DNA or RNA
strands as well as DNA or RNA duplexes or arrays of duplexes can be
found in /examples/USER/cgdna/util/. A technical report with more
found in examples/PACKAGES/cgdna/util/. A technical report with more
general information on the model, its implementation and performance
as well as the structure of the data and input file can be found
in /doc/src/PDF/USER-CGDNA.pdf.
in doc/src/PDF/CG-DNA.pdf.
IMPORTANT NOTE: This package can only be used if LAMMPS is compiled
with the MOLECULE and ASPHERE packages. These should be included in

View File

@ -0,0 +1,95 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "atom_vec_oxdna.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
AtomVecOxdna::AtomVecOxdna(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = Atom::MOLECULAR;
bonds_allow = 1;
mass_type = PER_TYPE;
atom->molecule_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = (char *) "id5p";
fields_copy = (char *) "id5p";
fields_comm = (char *) "";
fields_comm_vel = (char *) "";
fields_reverse = (char *) "";
fields_border = (char *) "id5p";
fields_border_vel = (char *) "";
fields_exchange = (char *) "id5p";
fields_restart = (char *) "id5p";
fields_create = (char *) "";
fields_data_atom = (char *) "id type x";
fields_data_vel = (char *) "id v";
setup_fields();
if(!force->newton_bond) error->warning(FLERR,"Write_data command requires newton on to preserve 3'->5' bond polarity");
}
/* ---------------------------------------------------------------------- */
AtomVecOxdna::~AtomVecOxdna() {}
/* ----------------------------------------------------------------------
set local copies of all grow ptrs used by this class, except defaults
needed in replicate when 2 atom classes exist and it calls pack_restart()
------------------------------------------------------------------------- */
void AtomVecOxdna::grow_pointers()
{
id5p = atom->id5p;
}
/* ----------------------------------------------------------------------
initialize atom quantity 5' partner
------------------------------------------------------------------------- */
void AtomVecOxdna::data_atom_post(int ilocal)
{
tagint *id5p = atom->id5p;
id5p[ilocal] = -1;
}
/* ----------------------------------------------------------------------
process bond information as per data file
store 5' partner to inform 3'->5' bond directionality
------------------------------------------------------------------------- */
void AtomVecOxdna::data_bonds_post(int /*m*/, int /*num_bond*/, tagint atom1, tagint atom2,
tagint id_offset)
{
int n;
tagint *id5p = atom->id5p;
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
}
if ((n = atom->map(atom1)) >= 0) { id5p[n] = atom2; }
}

View File

@ -1,4 +1,3 @@
// clang-format off
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -12,33 +11,37 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
#ifdef ATOM_CLASS
// clang-format off
FixStyle(reax/c/species/kk,FixReaxCSpeciesKokkos);
FixStyle(reax/c/species/kk/device,FixReaxCSpeciesKokkos);
FixStyle(reax/c/species/kk/host,FixReaxCSpeciesKokkos);
AtomStyle(oxdna,AtomVecOxdna);
// clang-format on
#else
#ifndef LMP_FIX_REAXC_SPECIES_KOKKOS_H
#define LMP_FIX_REAXC_SPECIES_KOKKOS_H
#ifndef LMP_ATOM_VEC_OXDNA_H
#define LMP_ATOM_VEC_OXDNA_H
#include "fix_reaxc_species.h"
#define BUFLEN 1000
#include "atom_vec.h"
namespace LAMMPS_NS {
class FixReaxCSpeciesKokkos : public FixReaxCSpecies {
class AtomVecOxdna : public AtomVec {
public:
FixReaxCSpeciesKokkos(class LAMMPS *, int, char **);
virtual ~FixReaxCSpeciesKokkos();
void init();
AtomVecOxdna(class LAMMPS *);
~AtomVecOxdna();
void grow_pointers();
void data_atom_post(int);
void data_bonds_post(int, int, tagint, tagint, tagint);
private:
void FindMolecule();
tagint *id5p;
};
}
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
*/

View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -19,30 +18,16 @@
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondOxdna2Fene::BondOxdna2Fene(LAMMPS *lmp) : BondOxdnaFene(lmp)
{
}
/* ---------------------------------------------------------------------- */
BondOxdna2Fene::~BondOxdna2Fene()
{
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxDNA2
------------------------------------------------------------------------- */
void BondOxdna2Fene::compute_interaction_sites(double e1[3], double e2[3],
double /*e3*/[3], double r[3])
void BondOxdna2Fene::compute_interaction_sites(double e1[3], double e2[3], double /*e3*/[3],
double r[3]) const
{
double d_cs_x=-0.34, d_cs_y=+0.3408;
r[0] = d_cs_x*e1[0] + d_cs_y*e2[0];
r[1] = d_cs_x*e1[1] + d_cs_y*e2[1];
r[2] = d_cs_x*e1[2] + d_cs_y*e2[2];
constexpr double d_cs_x = -0.34;
constexpr double d_cs_y = +0.3408;
r[0] = d_cs_x * e1[0] + d_cs_y * e2[0];
r[1] = d_cs_x * e1[1] + d_cs_y * e2[1];
r[2] = d_cs_x * e1[2] + d_cs_y * e2[2];
}

View File

@ -26,9 +26,9 @@ namespace LAMMPS_NS {
class BondOxdna2Fene : public BondOxdnaFene {
public:
BondOxdna2Fene(class LAMMPS *);
virtual ~BondOxdna2Fene();
virtual void compute_interaction_sites(double *, double *, double *, double *);
BondOxdna2Fene(class LAMMPS *lmp) : BondOxdnaFene(lmp) {}
virtual ~BondOxdna2Fene() {}
virtual void compute_interaction_sites(double *, double *, double *, double *) const;
};
} // namespace LAMMPS_NS

View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -17,7 +16,6 @@
#include "bond_oxdna_fene.h"
#include <cmath>
#include "atom.h"
#include "neighbor.h"
#include "comm.h"
@ -29,26 +27,19 @@
#include "atom_vec_ellipsoid.h"
#include "math_extra.h"
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondOxdnaFene::BondOxdnaFene(LAMMPS *lmp) : Bond(lmp)
{
}
/* ---------------------------------------------------------------------- */
BondOxdnaFene::~BondOxdnaFene()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(k);
memory->destroy(Delta);
memory->destroy(r0);
}
}
@ -57,14 +48,13 @@ BondOxdnaFene::~BondOxdnaFene()
compute vector COM-sugar-phosphate backbone interaction site in oxDNA
------------------------------------------------------------------------- */
void BondOxdnaFene::compute_interaction_sites(double e1[3], double /*e2*/[3],
double /*e3*/[3], double r[3])
double /*e3*/[3], double r[3]) const
{
double d_cs=-0.4;
constexpr double d_cs=-0.4;
r[0] = d_cs*e1[0];
r[1] = d_cs*e1[1];
r[2] = d_cs*e1[2];
}
/* ----------------------------------------------------------------------
@ -218,11 +208,8 @@ void BondOxdnaFene::compute(int eflag, int vflag)
// if r > 2*Delta something serious is wrong, abort
if (rlogarg < 0.1) {
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " %g",
update->ntimestep,atom->tag[a],atom->tag[b],r);
error->warning(FLERR,str);
error->warning(FLERR,"FENE bond too long: {} {} {} {}",
update->ntimestep,atom->tag[a],atom->tag[b],r);
rlogarg = 0.1;
}
@ -398,10 +385,8 @@ double BondOxdnaFene::single(int type, double rsq, int /*i*/, int /*j*/,
// if r > 2*Delta something serious is wrong, abort
if (rlogarg < 0.1) {
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %g",
update->ntimestep,sqrt(rsq));
error->warning(FLERR,str);
error->warning(FLERR,"FENE bond too long: {} {:.8}",
update->ntimestep,sqrt(rsq));
rlogarg = 0.1;
}

View File

@ -26,9 +26,9 @@ namespace LAMMPS_NS {
class BondOxdnaFene : public Bond {
public:
BondOxdnaFene(class LAMMPS *);
BondOxdnaFene(class LAMMPS *lmp) : Bond(lmp) {}
virtual ~BondOxdnaFene();
virtual void compute_interaction_sites(double *, double *, double *, double *);
virtual void compute_interaction_sites(double *, double *, double *, double *) const;
virtual void compute(int, int);
void coeff(int, char **);
void init_style();

View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -19,30 +18,17 @@
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondOxrna2Fene::BondOxrna2Fene(LAMMPS *lmp) : BondOxdnaFene(lmp)
{
}
/* ---------------------------------------------------------------------- */
BondOxrna2Fene::~BondOxrna2Fene()
{
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxRNA2
------------------------------------------------------------------------- */
void BondOxrna2Fene::compute_interaction_sites(double e1[3], double /*e2*/[3],
double e3[3], double r[3])
void BondOxrna2Fene::compute_interaction_sites(double e1[3], double /*e2*/[3], double e3[3],
double r[3]) const
{
double d_cs_x=-0.4, d_cs_z=+0.2;
r[0] = d_cs_x*e1[0] + d_cs_z*e3[0];
r[1] = d_cs_x*e1[1] + d_cs_z*e3[1];
r[2] = d_cs_x*e1[2] + d_cs_z*e3[2];
constexpr double d_cs_x = -0.4;
constexpr double d_cs_z = +0.2;
r[0] = d_cs_x * e1[0] + d_cs_z * e3[0];
r[1] = d_cs_x * e1[1] + d_cs_z * e3[1];
r[2] = d_cs_x * e1[2] + d_cs_z * e3[2];
}

View File

@ -26,9 +26,10 @@ namespace LAMMPS_NS {
class BondOxrna2Fene : public BondOxdnaFene {
public:
BondOxrna2Fene(class LAMMPS *);
virtual ~BondOxrna2Fene();
virtual void compute_interaction_sites(double *, double *, double *, double *);
BondOxrna2Fene(class LAMMPS *lmp) : BondOxdnaFene(lmp) {}
virtual ~BondOxrna2Fene() {}
virtual void compute_interaction_sites(double *, double *, double *, double *) const;
};
} // namespace LAMMPS_NS

View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -19,34 +18,19 @@
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairOxdna2Excv::PairOxdna2Excv(LAMMPS *lmp) : PairOxdnaExcv(lmp)
{
}
/* ---------------------------------------------------------------------- */
PairOxdna2Excv::~PairOxdna2Excv()
{
}
/* ----------------------------------------------------------------------
compute vector COM-excluded volume interaction sites in oxDNA2
------------------------------------------------------------------------- */
void PairOxdna2Excv::compute_interaction_sites(double e1[3], double e2[3],
double /*e3*/[3], double rs[3], double rb[3])
void PairOxdna2Excv::compute_interaction_sites(double e1[3], double e2[3], double /*e3*/[3],
double rs[3], double rb[3])
{
double d_cs_x=-0.34, d_cs_y=+0.3408, d_cb=+0.4;
double d_cs_x = -0.34, d_cs_y = +0.3408, d_cb = +0.4;
rs[0] = d_cs_x*e1[0] + d_cs_y*e2[0];
rs[1] = d_cs_x*e1[1] + d_cs_y*e2[1];
rs[2] = d_cs_x*e1[2] + d_cs_y*e2[2];
rb[0] = d_cb*e1[0];
rb[1] = d_cb*e1[1];
rb[2] = d_cb*e1[2];
rs[0] = d_cs_x * e1[0] + d_cs_y * e2[0];
rs[1] = d_cs_x * e1[1] + d_cs_y * e2[1];
rs[2] = d_cs_x * e1[2] + d_cs_y * e2[2];
rb[0] = d_cb * e1[0];
rb[1] = d_cb * e1[1];
rb[2] = d_cb * e1[2];
}

View File

@ -26,8 +26,8 @@ namespace LAMMPS_NS {
class PairOxdna2Excv : public PairOxdnaExcv {
public:
PairOxdna2Excv(class LAMMPS *);
virtual ~PairOxdna2Excv();
PairOxdna2Excv(class LAMMPS *lmp) : PairOxdnaExcv(lmp) {}
virtual ~PairOxdna2Excv() {}
virtual void compute_interaction_sites(double *, double *, double *, double *, double *);
};

View File

@ -17,19 +17,20 @@
#include "pair_oxdna_stk.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_oxdna.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_extra.h"
#include "memory.h"
#include "mf_oxdna.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
#include <utility>
#include "mf_oxdna.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "memory.h"
#include "error.h"
#include "atom_vec_ellipsoid.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace MFOxdna;
@ -42,26 +43,26 @@ PairOxdnaStk::PairOxdnaStk(LAMMPS *lmp) : Pair(lmp)
writedata = 1;
// sequence-specific stacking strength
// A:0 C:1 G:2 T:3, 5'- [i][j] -3'
// A:0 C:1 G:2 T:3, 3'- [i][j] -5'
eta_st[0][0] = 1.11960;
eta_st[0][1] = 1.00852;
eta_st[0][2] = 0.96950;
eta_st[0][3] = 0.99632;
eta_st[1][0] = 1.00852;
eta_st[2][0] = 0.96950;
eta_st[3][0] = 0.99632;
eta_st[1][0] = 1.01889;
eta_st[0][1] = 1.01889;
eta_st[1][1] = 0.97804;
eta_st[1][2] = 1.02681;
eta_st[1][3] = 0.96950;
eta_st[2][1] = 1.02681;
eta_st[3][1] = 0.96950;
eta_st[2][0] = 0.98169;
eta_st[2][1] = 1.05913;
eta_st[0][2] = 0.98169;
eta_st[1][2] = 1.05913;
eta_st[2][2] = 0.97804;
eta_st[2][3] = 1.00852;
eta_st[3][2] = 1.00852;
eta_st[3][0] = 0.94694;
eta_st[3][1] = 0.98169;
eta_st[3][2] = 1.01889;
eta_st[0][3] = 0.94694;
eta_st[1][3] = 0.98169;
eta_st[2][3] = 1.01889;
eta_st[3][3] = 0.96383;
}
@ -244,15 +245,16 @@ void PairOxdnaStk::compute(int eflag, int vflag)
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
tagint *id5p = atom->id5p;
AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int *ellipsoid = atom->ellipsoid;
int a,b,in,atype,btype;
int a,b,btemp,in,atype,btype;
double f1,f4t4,f4t5,f4t6,f5c1,f5c2;
double df1,df4t4,df4t5,df4t6,df5c1,df5c2;
double tptofp;
evdwl = 0.0;
ev_init(eflag,vflag);
@ -261,8 +263,19 @@ void PairOxdnaStk::compute(int eflag, int vflag)
for (in = 0; in < nbondlist; in++) {
a = bondlist[in][1];
b = bondlist[in][0];
a = bondlist[in][0];
b = bondlist[in][1];
// directionality test: a -> b is 3' -> 5'
if(atom->tag[b] != id5p[a]) {
btemp = b;
b = a;
a = btemp;
}
// a now in 3' direction, b in 5' direction
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
@ -279,30 +292,11 @@ void PairOxdnaStk::compute(int eflag, int vflag)
rb_cst[1] = d_cst*bx[1];
rb_cst[2] = d_cst*bx[2];
// vector stacking site b to a
delr_st[0] = x[a][0] + ra_cst[0] - x[b][0] - rb_cst[0];
delr_st[1] = x[a][1] + ra_cst[1] - x[b][1] - rb_cst[1];
delr_st[2] = x[a][2] + ra_cst[2] - x[b][2] - rb_cst[2];
// vector stacking site a to b
delr_st[0] = x[b][0] + rb_cst[0] - x[a][0] - ra_cst[0];
delr_st[1] = x[b][1] + rb_cst[1] - x[a][1] - ra_cst[1];
delr_st[2] = x[b][2] + rb_cst[2] - x[a][2] - ra_cst[2];
// test for directionality of vector b to a
tptofp = MFOxdna::is_3pto5p(delr_st,bz);
// if b to a is 5' to 3' we need to swap roles of a and b
if (tptofp == -1) {
std::swap(a,b);
std::swap(ax,bx);
std::swap(ay,by);
std::swap(az,bz);
std::swap(ra_cst,rb_cst);
delr_st[0] *= -1;
delr_st[1] *= -1;
delr_st[2] *= -1;
}
// a now in 5' direction, b in 3' direction
atype = type[a];
btype = type[b];
@ -325,9 +319,9 @@ void PairOxdnaStk::compute(int eflag, int vflag)
rb_cs[2] = d_cs*bx[2];
// vector backbone site b to a
delr_ss[0] = (x[a][0] + ra_cs[0] - x[b][0] - rb_cs[0]);
delr_ss[1] = (x[a][1] + ra_cs[1] - x[b][1] - rb_cs[1]);
delr_ss[2] = (x[a][2] + ra_cs[2] - x[b][2] - rb_cs[2]);
delr_ss[0] = (x[b][0] + rb_cs[0] - x[a][0] - ra_cs[0]);
delr_ss[1] = (x[b][1] + rb_cs[1] - x[a][1] - ra_cs[1]);
delr_ss[2] = (x[b][2] + rb_cs[2] - x[a][2] - ra_cs[2]);
rsq_ss = delr_ss[0]*delr_ss[0] + delr_ss[1]*delr_ss[1] + delr_ss[2]*delr_ss[2];
r_ss = sqrt(rsq_ss);
@ -345,7 +339,7 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (f1) {
// theta4 angle and correction
cost4 = MathExtra::dot3(az,bz);
cost4 = MathExtra::dot3(bz,az);
if (cost4 > 1.0) cost4 = 1.0;
if (cost4 < -1.0) cost4 = -1.0;
theta4 = acos(cost4);
@ -357,7 +351,7 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (f4t4) {
// theta5 angle and correction
cost5p = MathExtra::dot3(delr_st_norm,az);
cost5p = MathExtra::dot3(delr_st_norm,bz);
if (cost5p > 1.0) cost5p = 1.0;
if (cost5p < -1.0) cost5p = -1.0;
theta5p = acos(cost5p);
@ -368,16 +362,16 @@ void PairOxdnaStk::compute(int eflag, int vflag)
// early rejection criterium
if (f4t5) {
cost6p = MathExtra::dot3(delr_st_norm,bz);
cost6p = MathExtra::dot3(delr_st_norm,az);
if (cost6p > 1.0) cost6p = 1.0;
if (cost6p < -1.0) cost6p = -1.0;
theta6p = acos(cost6p);
cosphi1 = MathExtra::dot3(delr_ss_norm,ay);
cosphi1 = MathExtra::dot3(delr_ss_norm,by);
if (cosphi1 > 1.0) cosphi1 = 1.0;
if (cosphi1 < -1.0) cosphi1 = -1.0;
cosphi2 = MathExtra::dot3(delr_ss_norm,by);
cosphi2 = MathExtra::dot3(delr_ss_norm,ay);
if (cosphi2 > 1.0) cosphi2 = 1.0;
if (cosphi2 < -1.0) cosphi2 = -1.0;
@ -385,11 +379,10 @@ void PairOxdnaStk::compute(int eflag, int vflag)
b_st6[atype][btype], dtheta_st6_c[atype][btype]);
f5c1 = F5(-cosphi1, a_st1[atype][btype], -cosphi_st1_ast[atype][btype], b_st1[atype][btype],
cosphi_st1_c[atype][btype]);
-cosphi_st1_c[atype][btype]);
f5c2 = F5(-cosphi2, a_st2[atype][btype], -cosphi_st2_ast[atype][btype], b_st2[atype][btype],
cosphi_st2_c[atype][btype]);
-cosphi_st2_c[atype][btype]);
evdwl = f1 * f4t4 * f4t5 * f4t6 * f5c1 * f5c2;
@ -410,10 +403,10 @@ void PairOxdnaStk::compute(int eflag, int vflag)
b_st6[atype][btype], dtheta_st6_c[atype][btype])/sin(theta6p);
df5c1 = DF5(-cosphi1, a_st1[atype][btype], -cosphi_st1_ast[atype][btype], b_st1[atype][btype],
cosphi_st1_c[atype][btype]);
-cosphi_st1_c[atype][btype]);
df5c2 = DF5(-cosphi2, a_st2[atype][btype], -cosphi_st2_ast[atype][btype], b_st2[atype][btype],
cosphi_st2_c[atype][btype]);
-cosphi_st2_c[atype][btype]);
// force, torque and virial contribution for forces between stacking sites
@ -446,9 +439,9 @@ void PairOxdnaStk::compute(int eflag, int vflag)
finc = -f1 * f4t4 * df4t5 * f4t6 * f5c1 * f5c2 * rinv_st;
fpair += finc;
delf[0] += (delr_st_norm[0]*cost5p - az[0]) * finc;
delf[1] += (delr_st_norm[1]*cost5p - az[1]) * finc;
delf[2] += (delr_st_norm[2]*cost5p - az[2]) * finc;
delf[0] += (delr_st_norm[0]*cost5p - bz[0]) * finc;
delf[1] += (delr_st_norm[1]*cost5p - bz[1]) * finc;
delf[2] += (delr_st_norm[2]*cost5p - bz[2]) * finc;
}
@ -458,28 +451,27 @@ void PairOxdnaStk::compute(int eflag, int vflag)
finc = -f1 * f4t4 * f4t5 * df4t6 * f5c1 * f5c2 * rinv_st;
fpair += finc;
delf[0] += (delr_st_norm[0]*cost6p - bz[0]) * finc;
delf[1] += (delr_st_norm[1]*cost6p - bz[1]) * finc;
delf[2] += (delr_st_norm[2]*cost6p - bz[2]) * finc;
delf[0] += (delr_st_norm[0]*cost6p - az[0]) * finc;
delf[1] += (delr_st_norm[1]*cost6p - az[1]) * finc;
delf[2] += (delr_st_norm[2]*cost6p - az[2]) * finc;
}
// increment forces and torques
if (newton_bond || a < nlocal) {
f[a][0] += delf[0];
f[a][1] += delf[1];
f[a][2] += delf[2];
f[a][0] -= delf[0];
f[a][1] -= delf[1];
f[a][2] -= delf[2];
MathExtra::cross3(ra_cst,delf,delta);
}
if (newton_bond || b < nlocal) {
f[b][0] -= delf[0];
f[b][1] -= delf[1];
f[b][2] -= delf[2];
f[b][0] += delf[0];
f[b][1] += delf[1];
f[b][2] += delf[2];
MathExtra::cross3(rb_cst,delf,deltb);
@ -487,16 +479,16 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (newton_bond || a < nlocal) {
torque[a][0] += delta[0];
torque[a][1] += delta[1];
torque[a][2] += delta[2];
torque[a][0] -= delta[0];
torque[a][1] -= delta[1];
torque[a][2] -= delta[2];
}
if (newton_bond || b < nlocal) {
torque[b][0] -= deltb[0];
torque[b][1] -= deltb[1];
torque[b][2] -= deltb[2];
torque[b][0] += deltb[0];
torque[b][1] += deltb[1];
torque[b][2] += deltb[2];
}
@ -505,7 +497,7 @@ void PairOxdnaStk::compute(int eflag, int vflag)
// (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
if (evflag) ev_tally_xyz(a,b,nlocal,newton_bond,evdwl,
delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
delf[0],delf[1],delf[2],x[b][0]-x[a][0],x[b][1]-x[a][1],x[b][2]-x[a][2]);
// force, torque and virial contribution for forces between backbone sites
@ -529,9 +521,9 @@ void PairOxdnaStk::compute(int eflag, int vflag)
finc = -f1 * f4t4 * f4t5 * f4t6 * df5c1 * f5c2 * rinv_ss;
fpair += finc;
delf[0] += (delr_ss_norm[0]*cosphi1 - ay[0]) * finc;
delf[1] += (delr_ss_norm[1]*cosphi1 - ay[1]) * finc;
delf[2] += (delr_ss_norm[2]*cosphi1 - ay[2]) * finc;
delf[0] += (delr_ss_norm[0]*cosphi1 - by[0]) * finc;
delf[1] += (delr_ss_norm[1]*cosphi1 - by[1]) * finc;
delf[2] += (delr_ss_norm[2]*cosphi1 - by[2]) * finc;
}
@ -541,28 +533,27 @@ void PairOxdnaStk::compute(int eflag, int vflag)
finc = -f1 * f4t4 * f4t5 * f4t6 * f5c1 * df5c2 * rinv_ss;
fpair += finc;
delf[0] += (delr_ss_norm[0]*cosphi2 - by[0]) * finc;
delf[1] += (delr_ss_norm[1]*cosphi2 - by[1]) * finc;
delf[2] += (delr_ss_norm[2]*cosphi2 - by[2]) * finc;
delf[0] += (delr_ss_norm[0]*cosphi2 - ay[0]) * finc;
delf[1] += (delr_ss_norm[1]*cosphi2 - ay[1]) * finc;
delf[2] += (delr_ss_norm[2]*cosphi2 - ay[2]) * finc;
}
// increment forces and torques
if (newton_bond || a < nlocal) {
f[a][0] += delf[0];
f[a][1] += delf[1];
f[a][2] += delf[2];
f[a][0] -= delf[0];
f[a][1] -= delf[1];
f[a][2] -= delf[2];
MathExtra::cross3(ra_cs,delf,delta);
}
if (newton_bond || b < nlocal) {
f[b][0] -= delf[0];
f[b][1] -= delf[1];
f[b][2] -= delf[2];
f[b][0] += delf[0];
f[b][1] += delf[1];
f[b][2] += delf[2];
MathExtra::cross3(rb_cs,delf,deltb);
@ -570,22 +561,22 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (newton_bond || a < nlocal) {
torque[a][0] += delta[0];
torque[a][1] += delta[1];
torque[a][2] += delta[2];
torque[a][0] -= delta[0];
torque[a][1] -= delta[1];
torque[a][2] -= delta[2];
}
if (newton_bond || b < nlocal) {
torque[b][0] -= deltb[0];
torque[b][1] -= deltb[1];
torque[b][2] -= deltb[2];
torque[b][0] += deltb[0];
torque[b][1] += deltb[1];
torque[b][2] += deltb[2];
}
// increment virial only
if (evflag) ev_tally_xyz(a,b,nlocal,newton_bond,0.0,
delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
delf[0],delf[1],delf[2],x[b][0]-x[a][0],x[b][1]-x[a][1],x[b][2]-x[a][2]);
// pure torques not expressible as r x f
@ -600,7 +591,7 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (theta4) {
tpair = -f1 * df4t4 * f4t5 * f4t6 * f5c1 * f5c2;
MathExtra::cross3(bz,az,t4dir);
MathExtra::cross3(az,bz,t4dir);
delta[0] += t4dir[0]*tpair;
delta[1] += t4dir[1]*tpair;
@ -616,11 +607,11 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (theta5p) {
tpair = -f1 * f4t4 * df4t5 * f4t6 * f5c1 * f5c2;
MathExtra::cross3(delr_st_norm,az,t5pdir);
MathExtra::cross3(delr_st_norm,bz,t5pdir);
delta[0] += t5pdir[0] * tpair;
delta[1] += t5pdir[1] * tpair;
delta[2] += t5pdir[2] * tpair;
deltb[0] += t5pdir[0] * tpair;
deltb[1] += t5pdir[1] * tpair;
deltb[2] += t5pdir[2] * tpair;
}
@ -628,11 +619,11 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (theta6p) {
tpair = -f1 * f4t4 * f4t5 * df4t6 * f5c1 * f5c2;
MathExtra::cross3(delr_st_norm,bz,t6pdir);
MathExtra::cross3(delr_st_norm,az,t6pdir);
deltb[0] -= t6pdir[0] * tpair;
deltb[1] -= t6pdir[1] * tpair;
deltb[2] -= t6pdir[2] * tpair;
delta[0] -= t6pdir[0] * tpair;
delta[1] -= t6pdir[1] * tpair;
delta[2] -= t6pdir[2] * tpair;
}
@ -640,11 +631,11 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (cosphi1) {
tpair = -f1 * f4t4 * f4t5 * f4t6 * df5c1 * f5c2;
MathExtra::cross3(delr_ss_norm,ay,cosphi1dir);
MathExtra::cross3(delr_ss_norm,by,cosphi1dir);
delta[0] += cosphi1dir[0] * tpair;
delta[1] += cosphi1dir[1] * tpair;
delta[2] += cosphi1dir[2] * tpair;
deltb[0] += cosphi1dir[0] * tpair;
deltb[1] += cosphi1dir[1] * tpair;
deltb[2] += cosphi1dir[2] * tpair;
}
@ -652,27 +643,27 @@ void PairOxdnaStk::compute(int eflag, int vflag)
if (cosphi2) {
tpair = -f1 * f4t4 * f4t5 * f4t6 * f5c1 * df5c2;
MathExtra::cross3(delr_ss_norm,by,cosphi2dir);
MathExtra::cross3(delr_ss_norm,ay,cosphi2dir);
deltb[0] -= cosphi2dir[0] * tpair;
deltb[1] -= cosphi2dir[1] * tpair;
deltb[2] -= cosphi2dir[2] * tpair;
delta[0] -= cosphi2dir[0] * tpair;
delta[1] -= cosphi2dir[1] * tpair;
delta[2] -= cosphi2dir[2] * tpair;
}
// increment torques
if (newton_bond || a < nlocal) {
torque[a][0] += delta[0];
torque[a][1] += delta[1];
torque[a][2] += delta[2];
torque[a][0] -= delta[0];
torque[a][1] -= delta[1];
torque[a][2] -= delta[2];
}
if (newton_bond || b < nlocal) {
torque[b][0] -= deltb[0];
torque[b][1] -= deltb[1];
torque[b][2] -= deltb[2];
torque[b][0] += deltb[0];
torque[b][1] += deltb[1];
torque[b][2] += deltb[2];
}
@ -935,6 +926,17 @@ void PairOxdnaStk::coeff(int narg, char **arg)
}
/* ----------------------------------------------------------------------
atom_style hybrid bond ellipsoid oxdna required
------------------------------------------------------------------------- */
void PairOxdnaStk::init_style()
{
if (!atom->style_match("oxdna")) {
error->all(FLERR,"Must use 'atom_style hybrid bond ellipsoid oxdna' with pair style oxdna/stk, oxdna2/stk or oxrna2/stk");
}
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use regular
------------------------------------------------------------------------- */

View File

@ -32,6 +32,7 @@ class PairOxdnaStk : public Pair {
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);

View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -19,30 +18,15 @@
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairOxrna2Dh::PairOxrna2Dh(LAMMPS *lmp) : PairOxdna2Dh(lmp)
{
}
/* ---------------------------------------------------------------------- */
PairOxrna2Dh::~PairOxrna2Dh()
{
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxRNA2
------------------------------------------------------------------------- */
void PairOxrna2Dh::compute_interaction_sites(double e1[3], double /*e2*/[3],
double e3[3], double r[3])
void PairOxrna2Dh::compute_interaction_sites(double e1[3], double /*e2*/[3], double e3[3],
double r[3])
{
double d_cs_x=-0.4, d_cs_z=+0.2;
r[0] = d_cs_x*e1[0] + d_cs_z*e3[0];
r[1] = d_cs_x*e1[1] + d_cs_z*e3[1];
r[2] = d_cs_x*e1[2] + d_cs_z*e3[2];
double d_cs_x = -0.4, d_cs_z = +0.2;
r[0] = d_cs_x * e1[0] + d_cs_z * e3[0];
r[1] = d_cs_x * e1[1] + d_cs_z * e3[1];
r[2] = d_cs_x * e1[2] + d_cs_z * e3[2];
}

View File

@ -26,8 +26,8 @@ namespace LAMMPS_NS {
class PairOxrna2Dh : public PairOxdna2Dh {
public:
PairOxrna2Dh(class LAMMPS *);
virtual ~PairOxrna2Dh();
PairOxrna2Dh(class LAMMPS *lmp) : PairOxdna2Dh(lmp) {}
virtual ~PairOxrna2Dh() {}
virtual void compute_interaction_sites(double *, double *, double *, double *);
};

View File

@ -21,20 +21,6 @@
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairOxrna2Excv::PairOxrna2Excv(LAMMPS *lmp) : PairOxdnaExcv(lmp)
{
}
/* ---------------------------------------------------------------------- */
PairOxrna2Excv::~PairOxrna2Excv()
{
}
/* ----------------------------------------------------------------------
compute vector COM-excluded volume interaction sites in oxRNA2
------------------------------------------------------------------------- */

View File

@ -26,8 +26,8 @@ namespace LAMMPS_NS {
class PairOxrna2Excv : public PairOxdnaExcv {
public:
PairOxrna2Excv(class LAMMPS *);
virtual ~PairOxrna2Excv();
PairOxrna2Excv(class LAMMPS *lmp) : PairOxdnaExcv(lmp) {}
virtual ~PairOxrna2Excv() {}
virtual void compute_interaction_sites(double *, double *, double *, double *, double *);
};

View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -48,12 +47,4 @@ PairOxrna2Hbond::PairOxrna2Hbond(LAMMPS *lmp) : PairOxdnaHbond(lmp)
alpha_hb[3][1] = 1.00000;
alpha_hb[3][2] = 0.58655;
alpha_hb[3][3] = 1.00000;
}
/* ---------------------------------------------------------------------- */
PairOxrna2Hbond::~PairOxrna2Hbond()
{
}

View File

@ -27,7 +27,7 @@ namespace LAMMPS_NS {
class PairOxrna2Hbond : public PairOxdnaHbond {
public:
PairOxrna2Hbond(class LAMMPS *);
virtual ~PairOxrna2Hbond();
virtual ~PairOxrna2Hbond() {}
};
} // namespace LAMMPS_NS

View File

@ -19,6 +19,7 @@
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_oxdna.h"
#include "comm.h"
#include "error.h"
#include "force.h"
@ -44,26 +45,26 @@ PairOxrna2Stk::PairOxrna2Stk(LAMMPS *lmp) : Pair(lmp)
writedata = 1;
// sequence-specific stacking strength
// A:0 C:1 G:2 U:3, 5'- [i][j] -3'
// A:0 C:1 G:2 U:3, 3'- [i][j] -5'
eta_st[0][0] = 0.93851;
eta_st[0][1] = 1.12901;
eta_st[0][2] = 1.15626;
eta_st[0][3] = 0.88850;
eta_st[1][0] = 1.12901;
eta_st[2][0] = 1.15626;
eta_st[3][0] = 0.88850;
eta_st[1][0] = 0.86331;
eta_st[0][1] = 0.86331;
eta_st[1][1] = 1.05060;
eta_st[1][2] = 0.90982;
eta_st[1][3] = 0.83252;
eta_st[2][1] = 0.90982;
eta_st[3][1] = 0.83252;
eta_st[2][0] = 0.99407;
eta_st[2][1] = 1.14333;
eta_st[0][2] = 0.99407;
eta_st[1][2] = 1.14333;
eta_st[2][2] = 1.06573;
eta_st[2][3] = 0.91705;
eta_st[3][2] = 0.91705;
eta_st[3][0] = 0.98804;
eta_st[3][1] = 1.04949;
eta_st[3][2] = 1.12063;
eta_st[0][3] = 0.98804;
eta_st[1][3] = 1.04949;
eta_st[2][3] = 1.12063;
eta_st[3][3] = 0.83818;
}
@ -261,15 +262,16 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
tagint *id5p = atom->id5p;
AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int *ellipsoid = atom->ellipsoid;
int a,b,in,atype,btype;
int a,b,btemp,in,atype,btype;
double f1,f4t5,f4t6,f4t9,f4t10,f5c1,f5c2;
double df1,df4t5,df4t6,df4t9,df4t10,df5c1,df5c2;
double tptofp;
evdwl = 0.0;
ev_init(eflag,vflag);
@ -278,62 +280,43 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
for (in = 0; in < nbondlist; in++) {
a = bondlist[in][1];
b = bondlist[in][0];
a = bondlist[in][0];
b = bondlist[in][1];
// directionality test: a -> b is 3' -> 5'
if(atom->tag[b] != id5p[a]) {
btemp = b;
b = a;
a = btemp;
}
// a now in 3' direction, b in 5' direction
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
// vector COM a - 3'-stacking site a
ra_cst[0] = d_cst_x_3p*ax[0] + d_cst_y_3p*ay[0];
ra_cst[1] = d_cst_x_3p*ax[1] + d_cst_y_3p*ay[1];
ra_cst[2] = d_cst_x_3p*ax[2] + d_cst_y_3p*ay[2];
// vector COM a - 5'-stacking site a
ra_cst[0] = d_cst_x_5p*ax[0] + d_cst_y_5p*ay[0];
ra_cst[1] = d_cst_x_5p*ax[1] + d_cst_y_5p*ay[1];
ra_cst[2] = d_cst_x_5p*ax[2] + d_cst_y_5p*ay[2];
// vector COM b - 5'-stacking site b
rb_cst[0] = d_cst_x_5p*bx[0] + d_cst_y_5p*by[0];
rb_cst[1] = d_cst_x_5p*bx[1] + d_cst_y_5p*by[1];
rb_cst[2] = d_cst_x_5p*bx[2] + d_cst_y_5p*by[2];
// vector COM b - 3'-stacking site b
rb_cst[0] = d_cst_x_3p*bx[0] + d_cst_y_3p*by[0];
rb_cst[1] = d_cst_x_3p*bx[1] + d_cst_y_3p*by[1];
rb_cst[2] = d_cst_x_3p*bx[2] + d_cst_y_3p*by[2];
// vector 5'-stacking site b to 3'-stacking site a
delr_st[0] = x[a][0] + ra_cst[0] - x[b][0] - rb_cst[0];
delr_st[1] = x[a][1] + ra_cst[1] - x[b][1] - rb_cst[1];
delr_st[2] = x[a][2] + ra_cst[2] - x[b][2] - rb_cst[2];
// vector 5'-stacking site a to 3'-stacking site b
delr_st[0] = x[b][0] + rb_cst[0] - x[a][0] - ra_cst[0];
delr_st[1] = x[b][1] + rb_cst[1] - x[a][1] - ra_cst[1];
delr_st[2] = x[b][2] + rb_cst[2] - x[a][2] - ra_cst[2];
// test for directionality of vector b to a
tptofp = MFOxdna::is_3pto5p(delr_st,bz);
// if b to a is 5' to 3' we need to swap roles of a and b
if (tptofp == -1) {
std::swap(a,b);
std::swap(ax,bx);
std::swap(ay,by);
std::swap(az,bz);
}
// a now in 5' direction, b in 3' direction
atype = type[a];
btype = type[b];
// calculate again
// vector COM a - 3'-stacking site a
ra_cst[0] = d_cst_x_3p*ax[0] + d_cst_y_3p*ay[0];
ra_cst[1] = d_cst_x_3p*ax[1] + d_cst_y_3p*ay[1];
ra_cst[2] = d_cst_x_3p*ax[2] + d_cst_y_3p*ay[2];
// vector COM b - 5'-stacking site b
rb_cst[0] = d_cst_x_5p*bx[0] + d_cst_y_5p*by[0];
rb_cst[1] = d_cst_x_5p*bx[1] + d_cst_y_5p*by[1];
rb_cst[2] = d_cst_x_5p*bx[2] + d_cst_y_5p*by[2];
// vector 5'-stacking site b to 3'-stacking site a
delr_st[0] = x[a][0] + ra_cst[0] - x[b][0] - rb_cst[0];
delr_st[1] = x[a][1] + ra_cst[1] - x[b][1] - rb_cst[1];
delr_st[2] = x[a][2] + ra_cst[2] - x[b][2] - rb_cst[2];
rsq_st = delr_st[0]*delr_st[0] + delr_st[1]*delr_st[1] + delr_st[2]*delr_st[2];
r_st = sqrt(rsq_st);
rinv_st = 1.0/r_st;
@ -352,10 +335,10 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
rb_cs[1] = d_cs_x*bx[1] + d_cs_z*bz[1];
rb_cs[2] = d_cs_x*bx[2] + d_cs_z*bz[2];
// vector backbone site b to a
delr_ss[0] = (x[a][0] + ra_cs[0] - x[b][0] - rb_cs[0]);
delr_ss[1] = (x[a][1] + ra_cs[1] - x[b][1] - rb_cs[1]);
delr_ss[2] = (x[a][2] + ra_cs[2] - x[b][2] - rb_cs[2]);
// vector backbone site a to b
delr_ss[0] = (x[b][0] + rb_cs[0] - x[a][0] - ra_cs[0]);
delr_ss[1] = (x[b][1] + rb_cs[1] - x[a][1] - ra_cs[1]);
delr_ss[2] = (x[b][2] + rb_cs[2] - x[a][2] - ra_cs[2]);
rsq_ss = delr_ss[0]*delr_ss[0] + delr_ss[1]*delr_ss[1] + delr_ss[2]*delr_ss[2];
r_ss = sqrt(rsq_ss);
@ -373,7 +356,7 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
if (f1) {
// theta5 angle and correction
cost5p = MathExtra::dot3(delr_st_norm,az);
cost5p = MathExtra::dot3(delr_st_norm,bz);
if (cost5p > 1.0) cost5p = 1.0;
if (cost5p < -1.0) cost5p = -1.0;
theta5p = acos(cost5p);
@ -384,18 +367,18 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
// early rejection criterium
if (f4t5) {
cost6p = MathExtra::dot3(delr_st_norm,bz);
cost6p = MathExtra::dot3(delr_st_norm,az);
if (cost6p > 1.0) cost6p = 1.0;
if (cost6p < -1.0) cost6p = -1.0;
theta6p = acos(cost6p);
aux3p[0] = d3p_x * ax[0] + d3p_y * ay[0] + d3p_z * az[0];
aux3p[1] = d3p_x * ax[1] + d3p_y * ay[1] + d3p_z * az[1];
aux3p[2] = d3p_x * ax[2] + d3p_y * ay[2] + d3p_z * az[2];
aux5p[0] = d5p_x * ax[0] + d5p_y * ay[0] + d5p_z * az[0];
aux5p[1] = d5p_x * ax[1] + d5p_y * ay[1] + d5p_z * az[1];
aux5p[2] = d5p_x * ax[2] + d5p_y * ay[2] + d5p_z * az[2];
aux5p[0] = d5p_x * bx[0] + d5p_y * by[0] + d5p_z * bz[0];
aux5p[1] = d5p_x * bx[1] + d5p_y * by[1] + d5p_z * bz[1];
aux5p[2] = d5p_x * bx[2] + d5p_y * by[2] + d5p_z * bz[2];
aux3p[0] = d3p_x * bx[0] + d3p_y * by[0] + d3p_z * bz[0];
aux3p[1] = d3p_x * bx[1] + d3p_y * by[1] + d3p_z * bz[1];
aux3p[2] = d3p_x * bx[2] + d3p_y * by[2] + d3p_z * bz[2];
cost9 = MathExtra::dot3(delr_ss_norm,aux3p);
if (cost9 > 1.0) cost9 = 1.0;
@ -407,11 +390,11 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
if (cost10 < -1.0) cost10 = -1.0;
theta10 = acos(cost10);
cosphi1 = MathExtra::dot3(delr_ss_norm,ay);
cosphi1 = MathExtra::dot3(delr_ss_norm,by);
if (cosphi1 > 1.0) cosphi1 = 1.0;
if (cosphi1 < -1.0) cosphi1 = -1.0;
cosphi2 = MathExtra::dot3(delr_ss_norm,by);
cosphi2 = MathExtra::dot3(delr_ss_norm,ay);
if (cosphi2 > 1.0) cosphi2 = 1.0;
if (cosphi2 < -1.0) cosphi2 = -1.0;
@ -425,10 +408,10 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
b_st10[atype][btype], dtheta_st10_c[atype][btype]);
f5c1 = F5(-cosphi1, a_st1[atype][btype], -cosphi_st1_ast[atype][btype], b_st1[atype][btype],
cosphi_st1_c[atype][btype]);
-cosphi_st1_c[atype][btype]);
f5c2 = F5(-cosphi2, a_st2[atype][btype], -cosphi_st2_ast[atype][btype], b_st2[atype][btype],
cosphi_st2_c[atype][btype]);
-cosphi_st2_c[atype][btype]);
evdwl = f1 * f4t5 * f4t6 * f4t9 * f4t10 * f5c1 * f5c2;
@ -453,10 +436,10 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
b_st10[atype][btype], dtheta_st10_c[atype][btype])/sin(theta10);
df5c1 = DF5(-cosphi1, a_st1[atype][btype], -cosphi_st1_ast[atype][btype], b_st1[atype][btype],
cosphi_st1_c[atype][btype]);
-cosphi_st1_c[atype][btype]);
df5c2 = DF5(-cosphi2, a_st2[atype][btype], -cosphi_st2_ast[atype][btype], b_st2[atype][btype],
cosphi_st2_c[atype][btype]);
-cosphi_st2_c[atype][btype]);
// force, torque and virial contribution for forces between stacking sites
@ -489,9 +472,9 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
finc = -f1 * df4t5 * f4t6 * f4t9 * f4t10 * f5c1 * f5c2 * rinv_st;
fpair += finc;
delf[0] += (delr_st_norm[0]*cost5p - az[0]) * finc;
delf[1] += (delr_st_norm[1]*cost5p - az[1]) * finc;
delf[2] += (delr_st_norm[2]*cost5p - az[2]) * finc;
delf[0] += (delr_st_norm[0]*cost5p - bz[0]) * finc;
delf[1] += (delr_st_norm[1]*cost5p - bz[1]) * finc;
delf[2] += (delr_st_norm[2]*cost5p - bz[2]) * finc;
}
@ -501,28 +484,27 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
finc = -f1 * f4t5 * df4t6 * f4t9 * f4t10 * f5c1 * f5c2 * rinv_st;
fpair += finc;
delf[0] += (delr_st_norm[0]*cost6p - bz[0]) * finc;
delf[1] += (delr_st_norm[1]*cost6p - bz[1]) * finc;
delf[2] += (delr_st_norm[2]*cost6p - bz[2]) * finc;
delf[0] += (delr_st_norm[0]*cost6p - az[0]) * finc;
delf[1] += (delr_st_norm[1]*cost6p - az[1]) * finc;
delf[2] += (delr_st_norm[2]*cost6p - az[2]) * finc;
}
// increment forces and torques
if (newton_bond || a < nlocal) {
f[a][0] += delf[0];
f[a][1] += delf[1];
f[a][2] += delf[2];
f[a][0] -= delf[0];
f[a][1] -= delf[1];
f[a][2] -= delf[2];
MathExtra::cross3(ra_cst,delf,delta);
}
if (newton_bond || b < nlocal) {
f[b][0] -= delf[0];
f[b][1] -= delf[1];
f[b][2] -= delf[2];
f[b][0] += delf[0];
f[b][1] += delf[1];
f[b][2] += delf[2];
MathExtra::cross3(rb_cst,delf,deltb);
@ -530,16 +512,16 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
if (newton_bond || a < nlocal) {
torque[a][0] += delta[0];
torque[a][1] += delta[1];
torque[a][2] += delta[2];
torque[a][0] -= delta[0];
torque[a][1] -= delta[1];
torque[a][2] -= delta[2];
}
if (newton_bond || b < nlocal) {
torque[b][0] -= deltb[0];
torque[b][1] -= deltb[1];
torque[b][2] -= deltb[2];
torque[b][0] += deltb[0];
torque[b][1] += deltb[1];
torque[b][2] += deltb[2];
}
@ -547,8 +529,8 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
// NOTE: The virial is calculated on the 'molecular' basis.
// (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
if (evflag) ev_tally_xyz(a,b,nlocal,newton_bond,evdwl,
delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
if (evflag) ev_tally_xyz(b,a,nlocal,newton_bond,evdwl,
delf[0],delf[1],delf[2],x[b][0]-x[a][0],x[b][1]-x[a][1],x[b][2]-x[a][2]);
// force, torque and virial contribution for forces between backbone sites
@ -596,9 +578,9 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
finc = -f1 * f4t5 * f4t6 * f4t9 * f4t10 * df5c1 * f5c2 * rinv_ss;
fpair += finc;
delf[0] += (delr_ss_norm[0]*cosphi1 - ay[0]) * finc;
delf[1] += (delr_ss_norm[1]*cosphi1 - ay[1]) * finc;
delf[2] += (delr_ss_norm[2]*cosphi1 - ay[2]) * finc;
delf[0] += (delr_ss_norm[0]*cosphi1 - by[0]) * finc;
delf[1] += (delr_ss_norm[1]*cosphi1 - by[1]) * finc;
delf[2] += (delr_ss_norm[2]*cosphi1 - by[2]) * finc;
}
@ -608,28 +590,27 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
finc = -f1 * f4t5 * f4t6 * f4t9 * f4t10 * f5c1 * df5c2 * rinv_ss;
fpair += finc;
delf[0] += (delr_ss_norm[0]*cosphi2 - by[0]) * finc;
delf[1] += (delr_ss_norm[1]*cosphi2 - by[1]) * finc;
delf[2] += (delr_ss_norm[2]*cosphi2 - by[2]) * finc;
delf[0] += (delr_ss_norm[0]*cosphi2 - ay[0]) * finc;
delf[1] += (delr_ss_norm[1]*cosphi2 - ay[1]) * finc;
delf[2] += (delr_ss_norm[2]*cosphi2 - ay[2]) * finc;
}
// increment forces and torques
if (newton_bond || a < nlocal) {
f[a][0] += delf[0];
f[a][1] += delf[1];
f[a][2] += delf[2];
f[a][0] -= delf[0];
f[a][1] -= delf[1];
f[a][2] -= delf[2];
MathExtra::cross3(ra_cs,delf,delta);
}
if (newton_bond || b < nlocal) {
f[b][0] -= delf[0];
f[b][1] -= delf[1];
f[b][2] -= delf[2];
f[b][0] += delf[0];
f[b][1] += delf[1];
f[b][2] += delf[2];
MathExtra::cross3(rb_cs,delf,deltb);
@ -637,22 +618,22 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
if (newton_bond || a < nlocal) {
torque[a][0] += delta[0];
torque[a][1] += delta[1];
torque[a][2] += delta[2];
torque[a][0] -= delta[0];
torque[a][1] -= delta[1];
torque[a][2] -= delta[2];
}
if (newton_bond || b < nlocal) {
torque[b][0] -= deltb[0];
torque[b][1] -= deltb[1];
torque[b][2] -= deltb[2];
torque[b][0] += deltb[0];
torque[b][1] += deltb[1];
torque[b][2] += deltb[2];
}
// increment virial only
if (evflag) ev_tally_xyz(a,b,nlocal,newton_bond,0.0,
delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
if (evflag) ev_tally_xyz(b,a,nlocal,newton_bond,0.0,
delf[0],delf[1],delf[2],x[b][0]-x[a][0],x[b][1]-x[a][1],x[b][2]-x[a][2]);
// pure torques not expressible as r x f
@ -667,11 +648,11 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
if (theta5p) {
tpair = -f1 * df4t5 * f4t6 * f4t9 * f4t10 * f5c1 * f5c2;
MathExtra::cross3(delr_st_norm,az,t5pdir);
MathExtra::cross3(delr_st_norm,bz,t5pdir);
delta[0] += t5pdir[0] * tpair;
delta[1] += t5pdir[1] * tpair;
delta[2] += t5pdir[2] * tpair;
deltb[0] += t5pdir[0] * tpair;
deltb[1] += t5pdir[1] * tpair;
deltb[2] += t5pdir[2] * tpair;
}
@ -679,11 +660,11 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
if (theta6p) {
tpair = -f1 * f4t5 * df4t6 * f4t9 * f4t10 * f5c1 * f5c2;
MathExtra::cross3(delr_st_norm,bz,t6pdir);
MathExtra::cross3(delr_st_norm,az,t6pdir);
deltb[0] -= t6pdir[0] * tpair;
deltb[1] -= t6pdir[1] * tpair;
deltb[2] -= t6pdir[2] * tpair;
delta[0] -= t6pdir[0] * tpair;
delta[1] -= t6pdir[1] * tpair;
delta[2] -= t6pdir[2] * tpair;
}
@ -693,9 +674,9 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
tpair = -f1 * f4t5 * f4t6 * df4t9 * f4t10 * f5c1 * f5c2;
MathExtra::cross3(delr_ss_norm,aux3p,t9dir);
delta[0] += t9dir[0] * tpair;
delta[1] += t9dir[1] * tpair;
delta[2] += t9dir[2] * tpair;
deltb[0] += t9dir[0] * tpair;
deltb[1] += t9dir[1] * tpair;
deltb[2] += t9dir[2] * tpair;
}
@ -705,9 +686,9 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
tpair = -f1 * f4t5 * f4t6 * f4t9 * df4t10 * f5c1 * f5c2;
MathExtra::cross3(delr_ss_norm,aux5p,t10dir);
deltb[0] -= t10dir[0] * tpair;
deltb[1] -= t10dir[1] * tpair;
deltb[2] -= t10dir[2] * tpair;
delta[0] -= t10dir[0] * tpair;
delta[1] -= t10dir[1] * tpair;
delta[2] -= t10dir[2] * tpair;
}
@ -715,11 +696,11 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
if (cosphi1) {
tpair = -f1 * f4t5 * f4t6 * f4t9 * f4t10 * df5c1 * f5c2;
MathExtra::cross3(delr_ss_norm,ay,cosphi1dir);
MathExtra::cross3(delr_ss_norm,by,cosphi1dir);
delta[0] += cosphi1dir[0] * tpair;
delta[1] += cosphi1dir[1] * tpair;
delta[2] += cosphi1dir[2] * tpair;
deltb[0] += cosphi1dir[0] * tpair;
deltb[1] += cosphi1dir[1] * tpair;
deltb[2] += cosphi1dir[2] * tpair;
}
@ -727,27 +708,27 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
if (cosphi2) {
tpair = -f1 * f4t5 * f4t6 * f4t9 * f4t10 * f5c1 * df5c2;
MathExtra::cross3(delr_ss_norm,by,cosphi2dir);
MathExtra::cross3(delr_ss_norm,ay,cosphi2dir);
deltb[0] -= cosphi2dir[0] * tpair;
deltb[1] -= cosphi2dir[1] * tpair;
deltb[2] -= cosphi2dir[2] * tpair;
delta[0] -= cosphi2dir[0] * tpair;
delta[1] -= cosphi2dir[1] * tpair;
delta[2] -= cosphi2dir[2] * tpair;
}
// increment torques
if (newton_bond || a < nlocal) {
torque[a][0] += delta[0];
torque[a][1] += delta[1];
torque[a][2] += delta[2];
torque[a][0] -= delta[0];
torque[a][1] -= delta[1];
torque[a][2] -= delta[2];
}
if (newton_bond || b < nlocal) {
torque[b][0] -= deltb[0];
torque[b][1] -= deltb[1];
torque[b][2] -= deltb[2];
torque[b][0] += deltb[0];
torque[b][1] += deltb[1];
torque[b][2] += deltb[2];
}
@ -1027,6 +1008,17 @@ void PairOxrna2Stk::coeff(int narg, char **arg)
}
/* ----------------------------------------------------------------------
atom_style hybrid bond ellipsoid oxdna required
------------------------------------------------------------------------- */
void PairOxrna2Stk::init_style()
{
if (!atom->style_match("oxdna")) {
error->all(FLERR,"Must use 'atom_style hybrid bond ellipsoid oxdna' with pair style oxdna/stk, oxdna2/stk or oxrna2/stk");
}
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use regular
------------------------------------------------------------------------- */

View File

@ -31,6 +31,7 @@ class PairOxrna2Stk : public Pair {
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);

View File

@ -1,8 +1,9 @@
This package implements 3 commands which can be used in a LAMMPS input
This package implements 4 commands which can be used in a LAMMPS input
script:
pair_style lj/sdk
pair_style lj/sdk/coul/long
pair_style lj/sdk/coul/msm
angle_style sdk
These styles allow coarse grained MD simulations with the
@ -13,7 +14,7 @@ lipids and charged amino acids.
See the doc pages for these commands for details.
There are example scripts for using this package in
examples/USER/cgsdk
examples/PACKAGES/cgsdk
This is the second generation implementation reducing the the clutter
of the previous version. For many systems with long range

View File

@ -13,7 +13,7 @@ on building LAMMPS with external libraries. The settings in the
Makefile.lammps file in that directory must be correct for LAMMPS to
build correctly with this package installed.
The files in the USER-COLVARS package folder implement an interface
The files in the COLVARS package folder implement an interface
between LAMMPS and Colvars, originally written by Axel Kohlmeyer
(akohlmey@gmail.com) and maintained by Giacomo Fiorin
(giacomo.fiorin@gmail.com).
@ -36,7 +36,7 @@ which also includes citations to the articles documenting the various
methods that make use Colvars.
There are also example scripts for using this package in the folder
examples/USER/colvars, as well as the GitHub page for Colvars.
examples/PACKAGES/colvars, as well as the GitHub page for Colvars.
Please contact Giacomo Fiorin (giacomo.fiorin@gmail.com) for questions
regarding this package.

View File

@ -0,0 +1,3 @@
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2021-03-02"
#endif

View File

@ -480,7 +480,7 @@ void FixColvars::one_time_init()
memory->create(force_buf,3*num_coords,"colvars:force_buf");
if (me == 0) {
std::vector<int> &tl = *(proxy->modify_atom_ids());
std::vector<int> const &tl = *(proxy->get_atom_ids());
inthash_t *hashtable=new inthash_t;
inthash_init(hashtable, num_coords);
idmap = (void *)hashtable;
@ -563,7 +563,7 @@ void FixColvars::setup(int vflag)
if (me == 0) {
std::vector<int> &id = *(proxy->modify_atom_ids());
std::vector<int> const &id = *(proxy->get_atom_ids());
std::vector<int> &tp = *(proxy->modify_atom_types());
std::vector<cvm::atom_pos> &cd = *(proxy->modify_atom_positions());
std::vector<cvm::rvector> &of = *(proxy->modify_atom_total_forces());
@ -836,7 +836,7 @@ void FixColvars::post_force(int /*vflag*/)
if (me == 0) {
std::vector<cvm::rvector> &fo = *(proxy->modify_atom_new_colvar_forces());
std::vector<cvm::rvector> &fo = *(proxy->modify_atom_applied_forces());
double *fbuf = force_buf;
for (int j=0; j < num_coords; ++j) {

Some files were not shown because too many files have changed in this diff Show More