From d207710b43d4985774a35084ca0268c92b8b6127 Mon Sep 17 00:00:00 2001 From: Lewis M Russell <71833891+lrussell676@users.noreply.github.com> Date: Wed, 13 Oct 2021 10:47:41 +0100 Subject: [PATCH 01/13] lrefpos (nx/y/z) now working in oxdna_hbond. Being comm'ed and MPI working (#7) * use the term 'website' consistently (and not also 'web site') * update for clang-format * clarify * split off the programming/submission style guide to a separate file * Updates to support ROCm 4.3 in GPU package * update and reorder the description of the process for submitting contributions * correct and clarify Python compatibility * refactor style guide and integrate text from issue * update contribution guidelines for github * mention when testing may be added * integrate file with description of include file conventions * update github workflow doc * adapt section about domain decomposition from paper * use a more compact image * add communication section * update man page with missing flags and correct URLs * improve the load imbalance viz * add section about neighbor list construction * break large file into multiple smaller files by section and add toctree * fix typo * add section about parallelization in the OPENMP package * add -skipruin to help message * add discussion of OpenMP parallelization * spelling * add section on PPPM * use larger version of FFT grid comm image * minor tweak * Update compute angle doc page * Update Singularity definitions to use ROCm 4.3 * Update CUDA container definitions to CUDA 11.4 * Update container definitions to include PLUMED 2.7.2 * Update more definition files * RHEL8/CentOS8 PowerTools is now powertools * Add Rocky Linux 8 container definition * Add omega field to numpy_wrapper detection * Return None in case of null pointer * Add more atom fields in numpy_wrapper and correct csforce size * must not clear force array. will segfault in hybrid atom styles * update example for dynamically loading LAMMPS with current library API * update example to use current library interface. No need to use the namespace. * add note to README files about age of the example * simplify building shared libs on windows * detect a few more compilers * Revert "simplify building shared libs on windows" This reverts commit fa3429ab022b6bbb680011fee9a33f0a77839e1e. * step version strings for next patch release * fix mingw 32-bit vs 64-bit craziness * detect C++20 standard * build "fat" cuda binaries only with known toolkits * spelling * Try to improve the pair style hybrid docs This specifically tries to avoid the ambiguous use of "mixing" and clarify that similar is still different when pair styles are concerned. See discussion here: https://matsci.org/t/confusion-about-mixing-and-pair-coeff-section/38317/3 * spelling * use nullptr * use symbolic constant * small optimization * use cmath header instead of math.h * use explicit scoping when virtual dispatch is not available. * cosmetic changes * simplify/optimize code * simplify * Bugfix from Trung for crashes in pppm/gpu without local atoms * fix typo * refer to "XXX Coeffs" sections consistently * small tweaks from static code analysis * fix small bug * small tweaks * simplify and modernize code a little * use correct data type for MPI calls * simplify/modernize * remove dead code * about 1.5x speedup for pair style comb3 by using MathSpecial::powint() * small performance optimization for pair style comb * simplify * modernize * simplify * removed dead code, reformat * modernize * use explicit scoping when virtual dispatch is not (yet) available * reformat for increased readability * move misplaced #endif and make code more readable * make sure err_flag is initialized * modernize * remove redundant code: all struct members are initialized to defaults in the constructor * enforce initialization and thus silence compiler warnings * fix typo * provide more comprehensive suggestions for GPU neighbor list errors * add utils::logical() function to complement the *numeric() functions * Add stable link in docs * revert modernization change (for now) * remove unused variable * include EXTRA-DUMP in "most" * small tweaks * simplify * Add log file printing of KIM search directories in 'kim init' * use clang-format on kim_init.cpp * Improve style in response to Axel's suggestions * initialize all members * format changes * simplify. use utils::strdup() more. * small corrections * apply fix from balance command to fix balance * dead code removal * reformat strings * implement utils::current_date() convenience function to reduce replicated code * update list and order of include files from include-what-you-use analysis * handle changes in GAP repo * a few remaining updates to include statements * expand mapping to handle "style_*.h" header files correctly. * add support for compilation of OpenCL loader on FreeBSD * more iwyu header updates * small correction * fix typo * a few more (final?) IWYU updates * expand tests for numeric values * return int instead of bool to minimize code changes * fix spelling issues * some applications of the new function * fix typo * Change "offsite" to "external" to correct broken URLs to lammps.org * improve error message * insert missing atom-ID * convert yes/no on/off flags in the package command(s) * update version strings * update death tests for change in error message * correctly specify the destructor function name. * apply utils::logical() to more commands * apply utils::logical() in more places * for consistency with utils::logical() * only accept lower case to be consistent with the rest of the input * a few more converted commands and updates for unit tests * modernize and fix some memory leaks * adjust for compatibility with C++20 compilers * do not downgrade C++ standard when adding the KOKKOS package * undo "risky" C++20 related changes * mention how to set the path to the fftw3_omp library * correct paths to downloaded PACE package sources in lib * Update CMake variable descriptions * possible workaround for some GPU package neighbor list issue * final chunk of changes to apply utils::logical() * update suffix command unit tests * update citation info with new LAMMPS paper reference and acknowledge it * update some formulations as suggested by @sjplimp * add check and suitable error message when fp64 is required but not available * do not call memset on a null pointer * fix string formatting bugs in fix npt/cauchy * must use a soft core potential to avoid a singularity * in floating point math a*b may be zero even if both a>0 and b>0 * use proper integer type for atom IDs * Building voro++ lib as part of LAMMPS requires the "patch" program * silence output from hwloc when launching LAMMPS * detect double precision support according to OpenCL specs (1.2 and later) * Fix bug in Kokkos pair_eam_alloy * calling fwrite() with a null pointer causes undefined behavior. avoid it. * cosmetic * apply current include file conventions * include zstd libs in windows build * update .gitignore for recent additions * make check more obvious * step version strings for stable release * Adjust for kim-api bug * cleaner variant of version check, add directory numbering * hbond comm added for rsq_hb * rsq_hb removed, exyz added (no MPI comm yet) * Fix Colvars output files not written with "run 0" See: https://github.com/Colvars/colvars/commit/ff2f0d39ee5 which fixes a bug introduced in: https://github.com/Colvars/colvars/commit/1e964a542b The message applies to NAMD, but the logic used in LAMMPS when handling "run 0" is very similar. The Colvars version string is also updated, however this commit does not include other changes, such as the following: https://github.com/Colvars/colvars/pull/419 which were not fully completed before the LAMMPS Summer 2021 finalization. * add -std=c++11 to a number of machine makefiles for traditional make build * copy request to mention lammps.org form home page instructions for citing * be more specific about what the name of the LAMMPS executable can be also provide a few more examples without a machine suffix * small tweak * remove references to USER packages, have package lists alphabetically sorted "make package-update" or "make pu" must be processed in the special order because of inter-package dependencies * make "make package-update" and "make package-overwrite" less verbose * freeze versions of pip packages for processing the manual of the stable version this way we avoid surprises in case one of the packages get updated to an incompatible new version. these are know-to-work versions. * make sure the one_coeff flag is applied to sub-styles since the check for Pair::one_coeff was moved to the Input class (to reduce redundant code), hybrid substyles could "escape" that requirement. Thus checks have to be added to the hybrid coeff() methods. * Prevent neigh list from copying "unique" stencil/bin * compiling ML-HDNNP with downloaded n2p2 lib requires the sed command * detect and error out if BLAS/LAPACK libraries variables are a list This will cause external project compilation to fail since the semi-colons are converted to blanks, but one cannot properly escape the variables. So far the only viable solution seems to be to convert the scripts from using ExternalProject_add() to FetchContent and add_subdirectory() * portability improvement * must have patch command available to compile ScaFaCoS * only need Tcl not Tk to compile Tcl swig wrapper * correctly handle Tcl stub library if available * add missing keyword * hbond_pos added, MPI and values ok, Pair time slow. * make C library example work with strict C compilers * silence compiler warnings * make Nevery keyword per-reaction * recover cross-compilation with mingw64 * reverted wrong approach from last commits - now intpos flag - hbond_pos added - (a/b)xyz WiP * lrefpos working in serial, MPI wrong * attempt to merge doubles into n(xyz)[3] * save * Update pair_oxdna_hbond.cpp * hbond now working for MPI, comming lrefpos Co-authored-by: Axel Kohlmeyer Co-authored-by: Richard Berger Co-authored-by: Ryan S. Elliott Co-authored-by: Stan Gerald Moore Co-authored-by: Giacomo Fiorin Co-authored-by: Jacob Gissinger --- src/CG-DNA/bond_oxdna_fene.cpp | 3 +- src/CG-DNA/pair_oxdna_hbond.cpp | 178 ++++++++++++++++++++++++++++++-- src/CG-DNA/pair_oxdna_hbond.h | 8 ++ src/atom.cpp | 1 + src/atom.h | 4 + 5 files changed, 182 insertions(+), 12 deletions(-) diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp index 6b720f6a7f..6d92e4bcb0 100644 --- a/src/CG-DNA/bond_oxdna_fene.cpp +++ b/src/CG-DNA/bond_oxdna_fene.cpp @@ -255,7 +255,8 @@ void BondOxdnaFene::compute(int eflag, int vflag) if (evflag) ev_tally_xyz(a, b, nlocal, newton_bond, ebond, delf[0], delf[1], delf[2], x[a][0] - x[b][0], x[a][1] - x[b][1], x[a][2] - x[b][2]); - } + } + atom->lrefpos_flag = 0; // reset for next timestep } /* ---------------------------------------------------------------------- */ diff --git a/src/CG-DNA/pair_oxdna_hbond.cpp b/src/CG-DNA/pair_oxdna_hbond.cpp index 4bf0bc9c56..e268ffaded 100644 --- a/src/CG-DNA/pair_oxdna_hbond.cpp +++ b/src/CG-DNA/pair_oxdna_hbond.cpp @@ -62,6 +62,11 @@ PairOxdnaHbond::PairOxdnaHbond(LAMMPS *lmp) : Pair(lmp) alpha_hb[3][1] = 1.00000; alpha_hb[3][2] = 1.00000; alpha_hb[3][3] = 1.00000; + + // set comm size needed by this Pair + + comm_forward = 9; + comm_reverse = 9; } @@ -69,7 +74,12 @@ PairOxdnaHbond::PairOxdnaHbond(LAMMPS *lmp) : Pair(lmp) PairOxdnaHbond::~PairOxdnaHbond() { + if (allocated) { + + memory->destroy(nx); + memory->destroy(ny); + memory->destroy(nz); memory->destroy(setflag); memory->destroy(cutsq); @@ -149,9 +159,10 @@ void PairOxdnaHbond::compute(int eflag, int vflag) // vectors COM-h-bonding site in lab frame double ra_chb[3],rb_chb[3]; - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; + // Cartesian unit vectors in lab frame + // Only (a/b)x required here + double ax[3]; + double bx[3]; double **x = atom->x; double **f = atom->f; @@ -167,7 +178,7 @@ void PairOxdnaHbond::compute(int eflag, int vflag) AtomVecEllipsoid::Bonus *bonus = avec->bonus; int *ellipsoid = atom->ellipsoid; - int a,b,ia,ib,anum,bnum,atype,btype; + int n,a,b,in,ia,ib,anum,bnum,atype,btype; double f1,f4t1,f4t4,f4t2,f4t3,f4t7,f4t8; double df1,df4t1,df4t4,df4t2,df4t3,df4t7,df4t8; @@ -175,10 +186,41 @@ void PairOxdnaHbond::compute(int eflag, int vflag) evdwl = 0.0; ev_init(eflag,vflag); + nlocal = atom->nlocal; anum = list->inum; alist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; + + // loop over all local atoms, handle calculation of local reference frame + + if (atom->lrefpos_flag < atom->nmax) { + for (in = 0; in < atom->nlocal; in++) { + + int n = alist[in]; + + double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame + qn=bonus[ellipsoid[n]].quat; + MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp); + + nx[n][0] = nx_temp[0]; + nx[n][1] = nx_temp[1]; + nx[n][2] = nx_temp[2]; + ny[n][0] = ny_temp[0]; + ny[n][1] = ny_temp[1]; + ny[n][2] = ny_temp[2]; + nz[n][0] = nz_temp[0]; + nz[n][1] = nz_temp[1]; + nz[n][2] = nz_temp[2]; + + //printf("\n In top: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[n][0], nx[n][1], nx[n][2], atom->tag[n]); + + atom->lrefpos_flag += 1; + } + } + + //if (newton_pair) comm->reverse_comm_pair(this); + comm->forward_comm_pair(this); // loop over pair interaction neighbors of my atoms @@ -187,12 +229,17 @@ void PairOxdnaHbond::compute(int eflag, int vflag) a = alist[ia]; atype = type[a]; - qa=bonus[ellipsoid[a]].quat; - MathExtra::q_to_exyz(qa,ax,ay,az); + //printf("\n In A loop: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[a][0], nx[a][1], nx[a][2], atom->tag[a]); - ra_chb[0] = d_chb*ax[0]; + ax[0] = nx[a][0]; + ax[1] = nx[a][1]; + ax[2] = nx[a][2]; + + ra_chb[0] = d_chb*ax[0]; ra_chb[1] = d_chb*ax[1]; ra_chb[2] = d_chb*ax[2]; + + //printf("\n ax[0] = %f, ax[1] = %f, ax[2] = %f, id = %d", ax[0], ax[1], ax[2], atom->tag[a]); blist = firstneigh[a]; bnum = numneigh[a]; @@ -205,12 +252,15 @@ void PairOxdnaHbond::compute(int eflag, int vflag) btype = type[b]; - qb=bonus[ellipsoid[b]].quat; - MathExtra::q_to_exyz(qb,bx,by,bz); - - rb_chb[0] = d_chb*bx[0]; + bx[0] = nx[b][0]; + bx[1] = nx[b][1]; + bx[2] = nx[b][2]; + + rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; rb_chb[2] = d_chb*bx[2]; + + //printf("\n bx[0] = %f, bx[1] = %f, bx[2] = %f, id = %d", bx[0], bx[1], bx[2], atom->tag[b]); // vector h-bonding site b to a delr_hb[0] = x[a][0] + ra_chb[0] - x[b][0] - rb_chb[0]; @@ -264,6 +314,16 @@ void PairOxdnaHbond::compute(int eflag, int vflag) // early rejection criterium if (f4t3) { + + double az[3]; + double bz[3]; + + az[0] = nz[a][0]; + az[1] = nz[a][1]; + az[2] = nz[a][2]; + bz[0] = nz[b][0]; + bz[1] = nz[b][1]; + bz[2] = nz[b][2]; cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; @@ -560,6 +620,10 @@ void PairOxdnaHbond::allocate() for (int j = i; j <= n; j++) setflag[i][j] = 0; + memory->create(nx,atom->nmax,3,"pair:nx"); + memory->create(ny,atom->nmax,3,"pair:ny"); + memory->create(nz,atom->nmax,3,"pair:nz"); + memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(epsilon_hb,n+1,n+1,"pair:epsilon_hb"); @@ -1182,10 +1246,102 @@ void PairOxdnaHbond::write_data_all(FILE *fp) /* ---------------------------------------------------------------------- */ +int PairOxdnaHbond::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = nx[j][0]; + buf[m++] = nx[j][1]; + buf[m++] = nx[j][2]; + buf[m++] = ny[j][0]; + buf[m++] = ny[j][1]; + buf[m++] = ny[j][2]; + buf[m++] = nz[j][0]; + buf[m++] = nz[j][1]; + buf[m++] = nz[j][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void PairOxdnaHbond::unpack_forward_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + nx[i][0] = buf[m++]; + nx[i][1] = buf[m++]; + nx[i][2] = buf[m++]; + ny[i][0] = buf[m++]; + ny[i][1] = buf[m++]; + ny[i][2] = buf[m++]; + nz[i][0] = buf[m++]; + nz[i][1] = buf[m++]; + nz[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int PairOxdnaHbond::pack_reverse_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = nx[i][0]; + buf[m++] = nx[i][1]; + buf[m++] = nx[i][2]; + buf[m++] = ny[i][0]; + buf[m++] = ny[i][1]; + buf[m++] = ny[i][2]; + buf[m++] = nz[i][0]; + buf[m++] = nz[i][1]; + buf[m++] = nz[i][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void PairOxdnaHbond::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + nx[j][0] += buf[m++]; + nx[j][1] += buf[m++]; + nx[j][2] += buf[m++]; + ny[j][0] += buf[m++]; + ny[j][1] += buf[m++]; + ny[j][2] += buf[m++]; + nz[j][0] += buf[m++]; + nz[j][1] += buf[m++]; + nz[j][2] += buf[m++]; + + } +} + +/* ---------------------------------------------------------------------- */ + void *PairOxdnaHbond::extract(const char *str, int &dim) { dim = 2; + if (strcmp(str,"nx") == 0) return (void *) nx; + if (strcmp(str,"ny") == 0) return (void *) ny; + if (strcmp(str,"nz") == 0) return (void *) nz; + if (strcmp(str,"epsilon_hb") == 0) return (void *) epsilon_hb; if (strcmp(str,"a_hb") == 0) return (void *) a_hb; if (strcmp(str,"cut_hb_0") == 0) return (void *) cut_hb_0; diff --git a/src/CG-DNA/pair_oxdna_hbond.h b/src/CG-DNA/pair_oxdna_hbond.h index fdc8e9823b..e6dc47db56 100644 --- a/src/CG-DNA/pair_oxdna_hbond.h +++ b/src/CG-DNA/pair_oxdna_hbond.h @@ -41,6 +41,11 @@ class PairOxdnaHbond : public Pair { void write_data(FILE *); void write_data_all(FILE *); void *extract(const char *, int &); + + virtual int pack_forward_comm(int, int *, double *, int, int *); + virtual void unpack_forward_comm(int, int, double *); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); protected: // h-bonding interaction @@ -67,6 +72,9 @@ class PairOxdnaHbond : public Pair { double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast; double **b_hb8, **dtheta_hb8_c; + // per-atom arrays for q_to_exyz storage + double **nx, **ny, **nz; + int seqdepflag; virtual void allocate(); diff --git a/src/atom.cpp b/src/atom.cpp index 111ce7c93c..48adbf1eea 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -649,6 +649,7 @@ void Atom::set_atomflag_defaults() rho_flag = esph_flag = cv_flag = vest_flag = 0; dpd_flag = edpd_flag = tdpd_flag = 0; sp_flag = 0; + lrefpos_flag = 0; x0_flag = 0; smd_flag = damage_flag = 0; mesont_flag = 0; diff --git a/src/atom.h b/src/atom.h index e196d2d135..ee52c2c3d4 100644 --- a/src/atom.h +++ b/src/atom.h @@ -190,6 +190,10 @@ class Atom : protected Pointers { // SPIN package int sp_flag; + + // CG-DNA package + + int lrefpos_flag; // MACHDYN package From 20fec49635d78914ef4f6d5004529ac91376b2eb Mon Sep 17 00:00:00 2001 From: Lewis M Russell <71833891+lrussell676@users.noreply.github.com> Date: Tue, 26 Oct 2021 16:45:44 +0100 Subject: [PATCH 02/13] Intpos (#8) * use the term 'website' consistently (and not also 'web site') * update for clang-format * clarify * split off the programming/submission style guide to a separate file * Updates to support ROCm 4.3 in GPU package * update and reorder the description of the process for submitting contributions * correct and clarify Python compatibility * refactor style guide and integrate text from issue * update contribution guidelines for github * mention when testing may be added * integrate file with description of include file conventions * update github workflow doc * adapt section about domain decomposition from paper * use a more compact image * add communication section * update man page with missing flags and correct URLs * improve the load imbalance viz * add section about neighbor list construction * break large file into multiple smaller files by section and add toctree * fix typo * add section about parallelization in the OPENMP package * add -skipruin to help message * add discussion of OpenMP parallelization * spelling * add section on PPPM * use larger version of FFT grid comm image * minor tweak * Update compute angle doc page * Update Singularity definitions to use ROCm 4.3 * Update CUDA container definitions to CUDA 11.4 * Update container definitions to include PLUMED 2.7.2 * Update more definition files * RHEL8/CentOS8 PowerTools is now powertools * Add Rocky Linux 8 container definition * Add omega field to numpy_wrapper detection * Return None in case of null pointer * Add more atom fields in numpy_wrapper and correct csforce size * must not clear force array. will segfault in hybrid atom styles * update example for dynamically loading LAMMPS with current library API * update example to use current library interface. No need to use the namespace. * add note to README files about age of the example * simplify building shared libs on windows * detect a few more compilers * Revert "simplify building shared libs on windows" This reverts commit fa3429ab022b6bbb680011fee9a33f0a77839e1e. * step version strings for next patch release * fix mingw 32-bit vs 64-bit craziness * detect C++20 standard * build "fat" cuda binaries only with known toolkits * spelling * Try to improve the pair style hybrid docs This specifically tries to avoid the ambiguous use of "mixing" and clarify that similar is still different when pair styles are concerned. See discussion here: https://matsci.org/t/confusion-about-mixing-and-pair-coeff-section/38317/3 * spelling * use nullptr * use symbolic constant * small optimization * use cmath header instead of math.h * use explicit scoping when virtual dispatch is not available. * cosmetic changes * simplify/optimize code * simplify * Bugfix from Trung for crashes in pppm/gpu without local atoms * fix typo * refer to "XXX Coeffs" sections consistently * small tweaks from static code analysis * fix small bug * small tweaks * simplify and modernize code a little * use correct data type for MPI calls * simplify/modernize * remove dead code * about 1.5x speedup for pair style comb3 by using MathSpecial::powint() * small performance optimization for pair style comb * simplify * modernize * simplify * removed dead code, reformat * modernize * use explicit scoping when virtual dispatch is not (yet) available * reformat for increased readability * move misplaced #endif and make code more readable * make sure err_flag is initialized * modernize * remove redundant code: all struct members are initialized to defaults in the constructor * enforce initialization and thus silence compiler warnings * fix typo * provide more comprehensive suggestions for GPU neighbor list errors * add utils::logical() function to complement the *numeric() functions * Add stable link in docs * revert modernization change (for now) * remove unused variable * include EXTRA-DUMP in "most" * small tweaks * simplify * Add log file printing of KIM search directories in 'kim init' * use clang-format on kim_init.cpp * Improve style in response to Axel's suggestions * initialize all members * format changes * simplify. use utils::strdup() more. * small corrections * apply fix from balance command to fix balance * dead code removal * reformat strings * implement utils::current_date() convenience function to reduce replicated code * update list and order of include files from include-what-you-use analysis * handle changes in GAP repo * a few remaining updates to include statements * expand mapping to handle "style_*.h" header files correctly. * add support for compilation of OpenCL loader on FreeBSD * more iwyu header updates * small correction * fix typo * a few more (final?) IWYU updates * expand tests for numeric values * return int instead of bool to minimize code changes * fix spelling issues * some applications of the new function * fix typo * Change "offsite" to "external" to correct broken URLs to lammps.org * improve error message * insert missing atom-ID * convert yes/no on/off flags in the package command(s) * update version strings * update death tests for change in error message * correctly specify the destructor function name. * apply utils::logical() to more commands * apply utils::logical() in more places * for consistency with utils::logical() * only accept lower case to be consistent with the rest of the input * a few more converted commands and updates for unit tests * modernize and fix some memory leaks * adjust for compatibility with C++20 compilers * do not downgrade C++ standard when adding the KOKKOS package * undo "risky" C++20 related changes * mention how to set the path to the fftw3_omp library * correct paths to downloaded PACE package sources in lib * Update CMake variable descriptions * possible workaround for some GPU package neighbor list issue * final chunk of changes to apply utils::logical() * update suffix command unit tests * update citation info with new LAMMPS paper reference and acknowledge it * update some formulations as suggested by @sjplimp * add check and suitable error message when fp64 is required but not available * do not call memset on a null pointer * fix string formatting bugs in fix npt/cauchy * must use a soft core potential to avoid a singularity * in floating point math a*b may be zero even if both a>0 and b>0 * use proper integer type for atom IDs * Building voro++ lib as part of LAMMPS requires the "patch" program * silence output from hwloc when launching LAMMPS * detect double precision support according to OpenCL specs (1.2 and later) * Fix bug in Kokkos pair_eam_alloy * calling fwrite() with a null pointer causes undefined behavior. avoid it. * cosmetic * apply current include file conventions * include zstd libs in windows build * update .gitignore for recent additions * make check more obvious * step version strings for stable release * Adjust for kim-api bug * cleaner variant of version check, add directory numbering * hbond comm added for rsq_hb * rsq_hb removed, exyz added (no MPI comm yet) * Fix Colvars output files not written with "run 0" See: https://github.com/Colvars/colvars/commit/ff2f0d39ee5 which fixes a bug introduced in: https://github.com/Colvars/colvars/commit/1e964a542b The message applies to NAMD, but the logic used in LAMMPS when handling "run 0" is very similar. The Colvars version string is also updated, however this commit does not include other changes, such as the following: https://github.com/Colvars/colvars/pull/419 which were not fully completed before the LAMMPS Summer 2021 finalization. * add -std=c++11 to a number of machine makefiles for traditional make build * copy request to mention lammps.org form home page instructions for citing * be more specific about what the name of the LAMMPS executable can be also provide a few more examples without a machine suffix * small tweak * remove references to USER packages, have package lists alphabetically sorted "make package-update" or "make pu" must be processed in the special order because of inter-package dependencies * make "make package-update" and "make package-overwrite" less verbose * freeze versions of pip packages for processing the manual of the stable version this way we avoid surprises in case one of the packages get updated to an incompatible new version. these are know-to-work versions. * make sure the one_coeff flag is applied to sub-styles since the check for Pair::one_coeff was moved to the Input class (to reduce redundant code), hybrid substyles could "escape" that requirement. Thus checks have to be added to the hybrid coeff() methods. * Prevent neigh list from copying "unique" stencil/bin * compiling ML-HDNNP with downloaded n2p2 lib requires the sed command * detect and error out if BLAS/LAPACK libraries variables are a list This will cause external project compilation to fail since the semi-colons are converted to blanks, but one cannot properly escape the variables. So far the only viable solution seems to be to convert the scripts from using ExternalProject_add() to FetchContent and add_subdirectory() * portability improvement * must have patch command available to compile ScaFaCoS * only need Tcl not Tk to compile Tcl swig wrapper * correctly handle Tcl stub library if available * add missing keyword * hbond_pos added, MPI and values ok, Pair time slow. * make C library example work with strict C compilers * silence compiler warnings * make Nevery keyword per-reaction * recover cross-compilation with mingw64 * reverted wrong approach from last commits - now intpos flag - hbond_pos added - (a/b)xyz WiP * lrefpos working in serial, MPI wrong * attempt to merge doubles into n(xyz)[3] * save * Update pair_oxdna_hbond.cpp * hbond now working for MPI, comming lrefpos * extracting nxyz in excv/bond working Co-authored-by: Axel Kohlmeyer Co-authored-by: Richard Berger Co-authored-by: Ryan S. Elliott Co-authored-by: Stan Gerald Moore Co-authored-by: Giacomo Fiorin Co-authored-by: Jacob Gissinger Co-authored-by: Oliver Henrich --- src/CG-DNA/bond_oxdna_fene.cpp | 3 + src/CG-DNA/pair_oxdna_coaxstk.cpp | 2 + src/CG-DNA/pair_oxdna_excv.cpp | 124 ++++++++++++++++++++++++++++-- src/CG-DNA/pair_oxdna_excv.h | 6 ++ src/CG-DNA/pair_oxdna_hbond.cpp | 68 ++++++---------- src/CG-DNA/pair_oxdna_hbond.h | 10 +-- src/CG-DNA/pair_oxdna_stk.cpp | 2 + src/CG-DNA/pair_oxdna_xstk.cpp | 2 + 8 files changed, 155 insertions(+), 62 deletions(-) diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp index 6d92e4bcb0..468878c315 100644 --- a/src/CG-DNA/bond_oxdna_fene.cpp +++ b/src/CG-DNA/bond_oxdna_fene.cpp @@ -145,6 +145,9 @@ void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, doub ------------------------------------------------------------------------- */ void BondOxdnaFene::compute(int eflag, int vflag) { + + //printf("\n FENE-bond HERE, proc = %d \n", comm->me); + int a, b, in, type; double delf[3], delta[3], deltb[3]; // force, torque increment;; double delr[3], ebond, fbond; diff --git a/src/CG-DNA/pair_oxdna_coaxstk.cpp b/src/CG-DNA/pair_oxdna_coaxstk.cpp index 2715b7c5e2..c7ca9f9f47 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna_coaxstk.cpp @@ -106,6 +106,8 @@ PairOxdnaCoaxstk::~PairOxdnaCoaxstk() void PairOxdnaCoaxstk::compute(int eflag, int vflag) { + + //printf("\n Coaxstk HERE, proc = %d \n", comm->me); double delf[3],delt[3],delta[3],deltb[3]; // force, torque increment; double evdwl,fpair,finc,tpair,factor_lj; diff --git a/src/CG-DNA/pair_oxdna_excv.cpp b/src/CG-DNA/pair_oxdna_excv.cpp index d41976e86c..10eb371703 100644 --- a/src/CG-DNA/pair_oxdna_excv.cpp +++ b/src/CG-DNA/pair_oxdna_excv.cpp @@ -39,6 +39,9 @@ PairOxdnaExcv::PairOxdnaExcv(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; writedata = 1; + + // set comm size needed by this Pair + comm_forward = 9; } /* ---------------------------------------------------------------------- */ @@ -46,6 +49,10 @@ PairOxdnaExcv::PairOxdnaExcv(LAMMPS *lmp) : Pair(lmp) PairOxdnaExcv::~PairOxdnaExcv() { if (allocated) { + + memory->destroy(nx); + memory->destroy(ny); + memory->destroy(nz); memory->destroy(setflag); memory->destroy(cutsq); @@ -109,6 +116,8 @@ void PairOxdnaExcv::compute_interaction_sites(double e1[3], double /*e2*/[3], void PairOxdnaExcv::compute(int eflag, int vflag) { + //printf("\n ExcVol HERE, proc = %d \n", comm->me); + double delf[3],delta[3],deltb[3]; // force, torque increment; double evdwl,fpair,factor_lj; double rtmp_s[3],rtmp_b[3]; @@ -119,9 +128,8 @@ void PairOxdnaExcv::compute(int eflag, int vflag) double ra_cs[3],ra_cb[3]; double rb_cs[3],rb_cb[3]; - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double *special_lj = force->special_lj; double **x = atom->x; @@ -137,7 +145,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag) AtomVecEllipsoid::Bonus *bonus = avec->bonus; int *ellipsoid = atom->ellipsoid; - int a,b,ia,ib,anum,bnum,atype,btype; + int n,a,b,in,ia,ib,anum,bnum,atype,btype; evdwl = 0.0; ev_init(eflag,vflag); @@ -146,6 +154,36 @@ void PairOxdnaExcv::compute(int eflag, int vflag) alist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; + + // loop over all local atoms, handle calculation of local reference frame + + if (!atom->lrefpos_flag) { + for (in = 0; in < atom->nlocal; in++) { + + int n = alist[in]; + + double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame + qn=bonus[ellipsoid[n]].quat; + MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp); + + nx[n][0] = nx_temp[0]; + nx[n][1] = nx_temp[1]; + nx[n][2] = nx_temp[2]; + ny[n][0] = ny_temp[0]; + ny[n][1] = ny_temp[1]; + ny[n][2] = ny_temp[2]; + nz[n][0] = nz_temp[0]; + nz[n][1] = nz_temp[1]; + nz[n][2] = nz_temp[2]; + + //printf("\n In top: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[n][0], nx[n][1], nx[n][2], atom->tag[n]); + + atom->lrefpos_flag = 1; + } + } + + //if (newton_pair) comm->reverse_comm_pair(this); + comm->forward_comm_pair(this); // loop over pair interaction neighbors of my atoms @@ -154,8 +192,19 @@ void PairOxdnaExcv::compute(int eflag, int vflag) a = alist[ia]; atype = type[a]; - qa=bonus[ellipsoid[a]].quat; - MathExtra::q_to_exyz(qa,ax,ay,az); + //printf("\n In A loop: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[a][0], nx[a][1], nx[a][2], atom->tag[a]); + + ax[0] = nx[a][0]; + ax[1] = nx[a][1]; + ax[2] = nx[a][2]; + ay[0] = ny[a][0]; + ay[1] = ny[a][1]; + ay[2] = ny[a][2]; + az[0] = nz[a][0]; + az[1] = nz[a][1]; + az[2] = nz[a][2]; + + //printf("\n ax[0] = %f, ax[1] = %f, ax[2] = %f, id = %d", ax[0], ax[1], ax[2], atom->tag[a]); // vector COM - backbone and base site a compute_interaction_sites(ax,ay,az,ra_cs,ra_cb); @@ -179,8 +228,15 @@ void PairOxdnaExcv::compute(int eflag, int vflag) btype = type[b]; - qb=bonus[ellipsoid[b]].quat; - MathExtra::q_to_exyz(qb,bx,by,bz); + bx[0] = nx[b][0]; + bx[1] = nx[b][1]; + bx[2] = nx[b][2]; + by[0] = ny[b][0]; + by[1] = ny[b][1]; + by[2] = ny[b][2]; + bz[0] = nz[b][0]; + bz[1] = nz[b][1]; + bz[2] = nz[b][2]; // vector COM - backbone and base site b compute_interaction_sites(bx,by,bz,rb_cs,rb_cb); @@ -400,6 +456,10 @@ void PairOxdnaExcv::allocate() for (int j = i; j <= n; j++) setflag[i][j] = 0; + memory->create(nx,atom->nmax,3,"pair:nx"); + memory->create(ny,atom->nmax,3,"pair:ny"); + memory->create(nz,atom->nmax,3,"pair:nz"); + memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(epsilon_ss,n+1,n+1,"pair:epsilon_ss"); @@ -806,9 +866,57 @@ void PairOxdnaExcv::write_data_all(FILE *fp) /* ---------------------------------------------------------------------- */ +int PairOxdnaExcv::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = nx[j][0]; + buf[m++] = nx[j][1]; + buf[m++] = nx[j][2]; + buf[m++] = ny[j][0]; + buf[m++] = ny[j][1]; + buf[m++] = ny[j][2]; + buf[m++] = nz[j][0]; + buf[m++] = nz[j][1]; + buf[m++] = nz[j][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void PairOxdnaExcv::unpack_forward_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + nx[i][0] = buf[m++]; + nx[i][1] = buf[m++]; + nx[i][2] = buf[m++]; + ny[i][0] = buf[m++]; + ny[i][1] = buf[m++]; + ny[i][2] = buf[m++]; + nz[i][0] = buf[m++]; + nz[i][1] = buf[m++]; + nz[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + void *PairOxdnaExcv::extract(const char *str, int &dim) { dim = 2; + + if (strcmp(str,"nx") == 0) return (void *) nx; + if (strcmp(str,"ny") == 0) return (void *) ny; + if (strcmp(str,"nz") == 0) return (void *) nz; if (strcmp(str,"epsilon_ss") == 0) return (void *) epsilon_ss; if (strcmp(str,"sigma_ss") == 0) return (void *) sigma_ss; diff --git a/src/CG-DNA/pair_oxdna_excv.h b/src/CG-DNA/pair_oxdna_excv.h index c3ad6cf37a..1d2a871dd9 100644 --- a/src/CG-DNA/pair_oxdna_excv.h +++ b/src/CG-DNA/pair_oxdna_excv.h @@ -41,6 +41,9 @@ class PairOxdnaExcv : public Pair { void write_data(FILE *); void write_data_all(FILE *); void *extract(const char *, int &); + + virtual int pack_forward_comm(int, int *, double *, int, int *); + virtual void unpack_forward_comm(int, int, double *); protected: // s=sugar-phosphate backbone site, b=base site, st=stacking site @@ -52,6 +55,9 @@ class PairOxdnaExcv : public Pair { double **lj1_sb, **lj2_sb, **b_sb, **cut_sb_c, **cutsq_sb_c; double **epsilon_bb, **sigma_bb, **cut_bb_ast, **cutsq_bb_ast; double **lj1_bb, **lj2_bb, **b_bb, **cut_bb_c, **cutsq_bb_c; + + // per-atom arrays for q_to_exyz storage + double **nx, **ny, **nz; virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxdna_hbond.cpp b/src/CG-DNA/pair_oxdna_hbond.cpp index e268ffaded..083b60678e 100644 --- a/src/CG-DNA/pair_oxdna_hbond.cpp +++ b/src/CG-DNA/pair_oxdna_hbond.cpp @@ -192,35 +192,11 @@ void PairOxdnaHbond::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; - // loop over all local atoms, handle calculation of local reference frame - - if (atom->lrefpos_flag < atom->nmax) { - for (in = 0; in < atom->nlocal; in++) { - - int n = alist[in]; - - double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame - qn=bonus[ellipsoid[n]].quat; - MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp); - - nx[n][0] = nx_temp[0]; - nx[n][1] = nx_temp[1]; - nx[n][2] = nx_temp[2]; - ny[n][0] = ny_temp[0]; - ny[n][1] = ny_temp[1]; - ny[n][2] = ny_temp[2]; - nz[n][0] = nz_temp[0]; - nz[n][1] = nz_temp[1]; - nz[n][2] = nz_temp[2]; - - //printf("\n In top: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[n][0], nx[n][1], nx[n][2], atom->tag[n]); - - atom->lrefpos_flag += 1; - } - } - - //if (newton_pair) comm->reverse_comm_pair(this); - comm->forward_comm_pair(this); + // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + ny_xtrct = (double **) force->pair->extract("ny",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); // loop over pair interaction neighbors of my atoms @@ -231,11 +207,11 @@ void PairOxdnaHbond::compute(int eflag, int vflag) //printf("\n In A loop: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[a][0], nx[a][1], nx[a][2], atom->tag[a]); - ax[0] = nx[a][0]; - ax[1] = nx[a][1]; - ax[2] = nx[a][2]; + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; - ra_chb[0] = d_chb*ax[0]; + ra_chb[0] = d_chb*ax[0]; ra_chb[1] = d_chb*ax[1]; ra_chb[2] = d_chb*ax[2]; @@ -252,11 +228,11 @@ void PairOxdnaHbond::compute(int eflag, int vflag) btype = type[b]; - bx[0] = nx[b][0]; - bx[1] = nx[b][1]; - bx[2] = nx[b][2]; + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; - rb_chb[0] = d_chb*bx[0]; + rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; rb_chb[2] = d_chb*bx[2]; @@ -315,15 +291,15 @@ void PairOxdnaHbond::compute(int eflag, int vflag) // early rejection criterium if (f4t3) { - double az[3]; - double bz[3]; + double az[3]; + double bz[3]; - az[0] = nz[a][0]; - az[1] = nz[a][1]; - az[2] = nz[a][2]; - bz[0] = nz[b][0]; - bz[1] = nz[b][1]; - bz[2] = nz[b][2]; + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; @@ -1391,4 +1367,4 @@ void *PairOxdnaHbond::extract(const char *str, int &dim) if (strcmp(str,"dtheta_hb8_c") == 0) return (void *) dtheta_hb8_c; return nullptr; -} +} \ No newline at end of file diff --git a/src/CG-DNA/pair_oxdna_hbond.h b/src/CG-DNA/pair_oxdna_hbond.h index e6dc47db56..90ede130ba 100644 --- a/src/CG-DNA/pair_oxdna_hbond.h +++ b/src/CG-DNA/pair_oxdna_hbond.h @@ -41,11 +41,6 @@ class PairOxdnaHbond : public Pair { void write_data(FILE *); void write_data_all(FILE *); void *extract(const char *, int &); - - virtual int pack_forward_comm(int, int *, double *, int, int *); - virtual void unpack_forward_comm(int, int, double *); - int pack_reverse_comm(int, int, double *); - void unpack_reverse_comm(int, int *, double *); protected: // h-bonding interaction @@ -71,9 +66,8 @@ class PairOxdnaHbond : public Pair { double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast; double **b_hb8, **dtheta_hb8_c; - - // per-atom arrays for q_to_exyz storage - double **nx, **ny, **nz; + + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage int seqdepflag; diff --git a/src/CG-DNA/pair_oxdna_stk.cpp b/src/CG-DNA/pair_oxdna_stk.cpp index 2f1a0bf1b6..62de43c284 100644 --- a/src/CG-DNA/pair_oxdna_stk.cpp +++ b/src/CG-DNA/pair_oxdna_stk.cpp @@ -214,6 +214,8 @@ void PairOxdnaStk::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, void PairOxdnaStk::compute(int eflag, int vflag) { + + //printf("\n Stk HERE, proc = %d \n", comm->me); double delf[3],delta[3],deltb[3]; // force, torque increment; double evdwl,fpair,finc,tpair; diff --git a/src/CG-DNA/pair_oxdna_xstk.cpp b/src/CG-DNA/pair_oxdna_xstk.cpp index 447483e413..84e67666ee 100644 --- a/src/CG-DNA/pair_oxdna_xstk.cpp +++ b/src/CG-DNA/pair_oxdna_xstk.cpp @@ -112,6 +112,8 @@ PairOxdnaXstk::~PairOxdnaXstk() void PairOxdnaXstk::compute(int eflag, int vflag) { + //printf("\n Xstk HERE, proc = %d \n", comm->me); + double delf[3],delta[3],deltb[3]; // force, torque increment; double evdwl,fpair,finc,tpair,factor_lj; double delr_hb[3],delr_hb_norm[3],rsq_hb,r_hb,rinv_hb; From cf968ef286d24d07b977879231e0afd2cb7768ad Mon Sep 17 00:00:00 2001 From: Lewis M Russell <71833891+lrussell676@users.noreply.github.com> Date: Tue, 2 Nov 2021 09:52:56 +0000 Subject: [PATCH 03/13] Intpos (#10) * hbond comm added for rsq_hb * lrefpos removed, extract scaled for oxDNA1 * Update pair_oxdna_hbond.cpp --- src/CG-DNA/bond_oxdna_fene.cpp | 29 +++--- src/CG-DNA/bond_oxdna_fene.h | 1 + src/CG-DNA/pair_oxdna_coaxstk.cpp | 36 +++++-- src/CG-DNA/pair_oxdna_coaxstk.h | 2 + src/CG-DNA/pair_oxdna_excv.cpp | 53 +++++------ src/CG-DNA/pair_oxdna_hbond.cpp | 152 ++++-------------------------- src/CG-DNA/pair_oxdna_stk.cpp | 38 ++++++-- src/CG-DNA/pair_oxdna_stk.h | 3 + src/CG-DNA/pair_oxdna_xstk.cpp | 32 +++++-- src/CG-DNA/pair_oxdna_xstk.h | 3 + src/atom.cpp | 1 - src/atom.h | 4 - 12 files changed, 148 insertions(+), 206 deletions(-) diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp index 468878c315..c46556768c 100644 --- a/src/CG-DNA/bond_oxdna_fene.cpp +++ b/src/CG-DNA/bond_oxdna_fene.cpp @@ -26,6 +26,7 @@ #include "atom_vec_ellipsoid.h" #include "math_extra.h" +#include "pair.h" #include @@ -144,9 +145,7 @@ void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, doub s=sugar-phosphate backbone site, b=base site, st=stacking site ------------------------------------------------------------------------- */ void BondOxdnaFene::compute(int eflag, int vflag) -{ - - //printf("\n FENE-bond HERE, proc = %d \n", comm->me); +{ int a, b, in, type; double delf[3], delta[3], deltb[3]; // force, torque increment;; @@ -155,9 +154,9 @@ void BondOxdnaFene::compute(int eflag, int vflag) double r, rr0, rr0sq; // vectors COM-backbone site in lab frame double ra_cs[3], rb_cs[3]; - - double *qa, ax[3], ay[3], az[3]; - double *qb, bx[3], by[3], bz[3]; + // Cartesian unit vectors in lab frame + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -174,6 +173,12 @@ void BondOxdnaFene::compute(int eflag, int vflag) ebond = 0.0; ev_init(eflag, vflag); + + // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + ny_xtrct = (double **) force->pair->extract("ny",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); // loop over FENE bonds @@ -183,10 +188,13 @@ void BondOxdnaFene::compute(int eflag, int vflag) b = bondlist[in][0]; type = bondlist[in][2]; - 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); + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + //(a/b)y/z not needed here as oxDNA(1) co-linear // vector COM-backbone site a and b compute_interaction_sites(ax, ay, az, ra_cs); @@ -259,7 +267,6 @@ void BondOxdnaFene::compute(int eflag, int vflag) ev_tally_xyz(a, b, nlocal, newton_bond, ebond, delf[0], delf[1], delf[2], x[a][0] - x[b][0], x[a][1] - x[b][1], x[a][2] - x[b][2]); } - atom->lrefpos_flag = 0; // reset for next timestep } /* ---------------------------------------------------------------------- */ diff --git a/src/CG-DNA/bond_oxdna_fene.h b/src/CG-DNA/bond_oxdna_fene.h index 42b542a6fb..b24757781f 100644 --- a/src/CG-DNA/bond_oxdna_fene.h +++ b/src/CG-DNA/bond_oxdna_fene.h @@ -40,6 +40,7 @@ class BondOxdnaFene : public Bond { protected: double *k, *Delta, *r0; // FENE + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage void allocate(); void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double); diff --git a/src/CG-DNA/pair_oxdna_coaxstk.cpp b/src/CG-DNA/pair_oxdna_coaxstk.cpp index c7ca9f9f47..f76ad33147 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna_coaxstk.cpp @@ -131,10 +131,9 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) // vectors COM-backbone site, COM-stacking site in lab frame double ra_cs[3],ra_cst[3]; double rb_cs[3],rb_cst[3]; - - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; + // Cartesian unit vectors in lab frame + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -162,6 +161,12 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) alist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; + + // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + ny_xtrct = (double **) force->pair->extract("ny",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); // loop over pair interaction neighbors of my atoms @@ -170,8 +175,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) a = alist[ia]; atype = type[a]; - qa=bonus[ellipsoid[a]].quat; - MathExtra::q_to_exyz(qa,ax,ay,az); + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + //a(y/z) not needed here as oxDNA(1) co-linear // vector COM a - stacking site a ra_cst[0] = d_cst*ax[0]; @@ -194,8 +201,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) btype = type[b]; - qb=bonus[ellipsoid[b]].quat; - MathExtra::q_to_exyz(qb,bx,by,bz); + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + //b(y/z) not needed here as oxDNA(1) co-linear // vector COM b - stacking site b rb_cst[0] = d_cst*bx[0]; @@ -246,6 +255,13 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) // early rejection criterium if (f4t1) { + + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; @@ -382,6 +398,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) // cosphi3 and cosphi4 (=cosphi3) force and virial if (cosphi3) { + ay[0] = ny_xtrct[a][0]; + ay[1] = ny_xtrct[a][1]; + ay[2] = ny_xtrct[a][2]; + finc = -f2 * f4t1* f4t4 * f4t5 * f4t6 * 2.0 * f5c3 * df5c3 * factor_lj; fpair += finc; diff --git a/src/CG-DNA/pair_oxdna_coaxstk.h b/src/CG-DNA/pair_oxdna_coaxstk.h index 5771bbe592..d634ccdf00 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.h +++ b/src/CG-DNA/pair_oxdna_coaxstk.h @@ -62,6 +62,8 @@ class PairOxdnaCoaxstk : public Pair { double **a_cxst3p, **cosphi_cxst3p_ast, **b_cxst3p, **cosphi_cxst3p_c; double **a_cxst4p, **cosphi_cxst4p_ast, **b_cxst4p, **cosphi_cxst4p_c; + + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxdna_excv.cpp b/src/CG-DNA/pair_oxdna_excv.cpp index 10eb371703..2edad16c83 100644 --- a/src/CG-DNA/pair_oxdna_excv.cpp +++ b/src/CG-DNA/pair_oxdna_excv.cpp @@ -127,9 +127,10 @@ void PairOxdnaExcv::compute(int eflag, int vflag) // vectors COM-backbone site, COM-base site in lab frame double ra_cs[3],ra_cb[3]; double rb_cs[3],rb_cb[3]; - + // Cartesian unit vectors in lab frame double ax[3],ay[3],az[3]; double bx[3],by[3],bz[3]; + double *special_lj = force->special_lj; double **x = atom->x; @@ -157,32 +158,26 @@ void PairOxdnaExcv::compute(int eflag, int vflag) // loop over all local atoms, handle calculation of local reference frame - if (!atom->lrefpos_flag) { - for (in = 0; in < atom->nlocal; in++) { - - int n = alist[in]; - - double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame - qn=bonus[ellipsoid[n]].quat; - MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp); - - nx[n][0] = nx_temp[0]; - nx[n][1] = nx_temp[1]; - nx[n][2] = nx_temp[2]; - ny[n][0] = ny_temp[0]; - ny[n][1] = ny_temp[1]; - ny[n][2] = ny_temp[2]; - nz[n][0] = nz_temp[0]; - nz[n][1] = nz_temp[1]; - nz[n][2] = nz_temp[2]; - - //printf("\n In top: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[n][0], nx[n][1], nx[n][2], atom->tag[n]); - - atom->lrefpos_flag = 1; - } + for (in = 0; in < atom->nlocal; in++) { + + int n = alist[in]; + + double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame + qn=bonus[ellipsoid[n]].quat; + MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp); + + nx[n][0] = nx_temp[0]; + nx[n][1] = nx_temp[1]; + nx[n][2] = nx_temp[2]; + ny[n][0] = ny_temp[0]; + ny[n][1] = ny_temp[1]; + ny[n][2] = ny_temp[2]; + nz[n][0] = nz_temp[0]; + nz[n][1] = nz_temp[1]; + nz[n][2] = nz_temp[2]; + } - //if (newton_pair) comm->reverse_comm_pair(this); comm->forward_comm_pair(this); // loop over pair interaction neighbors of my atoms @@ -190,9 +185,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag) for (ia = 0; ia < anum; ia++) { a = alist[ia]; - atype = type[a]; - - //printf("\n In A loop: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[a][0], nx[a][1], nx[a][2], atom->tag[a]); + atype = type[a]; ax[0] = nx[a][0]; ax[1] = nx[a][1]; @@ -202,9 +195,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag) ay[2] = ny[a][2]; az[0] = nz[a][0]; az[1] = nz[a][1]; - az[2] = nz[a][2]; - - //printf("\n ax[0] = %f, ax[1] = %f, ax[2] = %f, id = %d", ax[0], ax[1], ax[2], atom->tag[a]); + az[2] = nz[a][2]; // vector COM - backbone and base site a compute_interaction_sites(ax,ay,az,ra_cs,ra_cb); diff --git a/src/CG-DNA/pair_oxdna_hbond.cpp b/src/CG-DNA/pair_oxdna_hbond.cpp index 083b60678e..f4baa74787 100644 --- a/src/CG-DNA/pair_oxdna_hbond.cpp +++ b/src/CG-DNA/pair_oxdna_hbond.cpp @@ -62,11 +62,6 @@ PairOxdnaHbond::PairOxdnaHbond(LAMMPS *lmp) : Pair(lmp) alpha_hb[3][1] = 1.00000; alpha_hb[3][2] = 1.00000; alpha_hb[3][3] = 1.00000; - - // set comm size needed by this Pair - - comm_forward = 9; - comm_reverse = 9; } @@ -74,12 +69,7 @@ PairOxdnaHbond::PairOxdnaHbond(LAMMPS *lmp) : Pair(lmp) PairOxdnaHbond::~PairOxdnaHbond() { - if (allocated) { - - memory->destroy(nx); - memory->destroy(ny); - memory->destroy(nz); memory->destroy(setflag); memory->destroy(cutsq); @@ -158,11 +148,9 @@ void PairOxdnaHbond::compute(int eflag, int vflag) double d_chb=+0.4; // vectors COM-h-bonding site in lab frame double ra_chb[3],rb_chb[3]; - // Cartesian unit vectors in lab frame - // Only (a/b)x required here - double ax[3]; - double bx[3]; + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -178,7 +166,7 @@ void PairOxdnaHbond::compute(int eflag, int vflag) AtomVecEllipsoid::Bonus *bonus = avec->bonus; int *ellipsoid = atom->ellipsoid; - int n,a,b,in,ia,ib,anum,bnum,atype,btype; + int a,b,ia,ib,anum,bnum,atype,btype; double f1,f4t1,f4t4,f4t2,f4t3,f4t7,f4t8; double df1,df4t1,df4t4,df4t2,df4t3,df4t7,df4t8; @@ -186,7 +174,6 @@ void PairOxdnaHbond::compute(int eflag, int vflag) evdwl = 0.0; ev_init(eflag,vflag); - nlocal = atom->nlocal; anum = list->inum; alist = list->ilist; numneigh = list->numneigh; @@ -205,17 +192,13 @@ void PairOxdnaHbond::compute(int eflag, int vflag) a = alist[ia]; atype = type[a]; - //printf("\n In A loop: nx[0] = %f, nx[1] = %f, nx[2] = %f, id = %d", nx[a][0], nx[a][1], nx[a][2], atom->tag[a]); - ax[0] = nx_xtrct[a][0]; - ax[1] = nx_xtrct[a][1]; - ax[2] = nx_xtrct[a][2]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; - ra_chb[0] = d_chb*ax[0]; + ra_chb[0] = d_chb*ax[0]; ra_chb[1] = d_chb*ax[1]; ra_chb[2] = d_chb*ax[2]; - - //printf("\n ax[0] = %f, ax[1] = %f, ax[2] = %f, id = %d", ax[0], ax[1], ax[2], atom->tag[a]); blist = firstneigh[a]; bnum = numneigh[a]; @@ -228,16 +211,14 @@ void PairOxdnaHbond::compute(int eflag, int vflag) btype = type[b]; - bx[0] = nx_xtrct[b][0]; - bx[1] = nx_xtrct[b][1]; - bx[2] = nx_xtrct[b][2]; + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; - rb_chb[0] = d_chb*bx[0]; + rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; rb_chb[2] = d_chb*bx[2]; - //printf("\n bx[0] = %f, bx[1] = %f, bx[2] = %f, id = %d", bx[0], bx[1], bx[2], atom->tag[b]); - // vector h-bonding site b to a delr_hb[0] = x[a][0] + ra_chb[0] - x[b][0] - rb_chb[0]; delr_hb[1] = x[a][1] + ra_chb[1] - x[b][1] - rb_chb[1]; @@ -289,17 +270,14 @@ void PairOxdnaHbond::compute(int eflag, int vflag) b_hb3[atype][btype], dtheta_hb3_c[atype][btype]); // early rejection criterium - if (f4t3) { + if (f4t3) { - double az[3]; - double bz[3]; - - az[0] = nz_xtrct[a][0]; - az[1] = nz_xtrct[a][1]; - az[2] = nz_xtrct[a][2]; - bz[0] = nz_xtrct[b][0]; - bz[1] = nz_xtrct[b][1]; - bz[2] = nz_xtrct[b][2]; + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; @@ -596,10 +574,6 @@ void PairOxdnaHbond::allocate() for (int j = i; j <= n; j++) setflag[i][j] = 0; - memory->create(nx,atom->nmax,3,"pair:nx"); - memory->create(ny,atom->nmax,3,"pair:ny"); - memory->create(nz,atom->nmax,3,"pair:nz"); - memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(epsilon_hb,n+1,n+1,"pair:epsilon_hb"); @@ -1222,102 +1196,10 @@ void PairOxdnaHbond::write_data_all(FILE *fp) /* ---------------------------------------------------------------------- */ -int PairOxdnaHbond::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = nx[j][0]; - buf[m++] = nx[j][1]; - buf[m++] = nx[j][2]; - buf[m++] = ny[j][0]; - buf[m++] = ny[j][1]; - buf[m++] = ny[j][2]; - buf[m++] = nz[j][0]; - buf[m++] = nz[j][1]; - buf[m++] = nz[j][2]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void PairOxdnaHbond::unpack_forward_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - nx[i][0] = buf[m++]; - nx[i][1] = buf[m++]; - nx[i][2] = buf[m++]; - ny[i][0] = buf[m++]; - ny[i][1] = buf[m++]; - ny[i][2] = buf[m++]; - nz[i][0] = buf[m++]; - nz[i][1] = buf[m++]; - nz[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int PairOxdnaHbond::pack_reverse_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = nx[i][0]; - buf[m++] = nx[i][1]; - buf[m++] = nx[i][2]; - buf[m++] = ny[i][0]; - buf[m++] = ny[i][1]; - buf[m++] = ny[i][2]; - buf[m++] = nz[i][0]; - buf[m++] = nz[i][1]; - buf[m++] = nz[i][2]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void PairOxdnaHbond::unpack_reverse_comm(int n, int *list, double *buf) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - nx[j][0] += buf[m++]; - nx[j][1] += buf[m++]; - nx[j][2] += buf[m++]; - ny[j][0] += buf[m++]; - ny[j][1] += buf[m++]; - ny[j][2] += buf[m++]; - nz[j][0] += buf[m++]; - nz[j][1] += buf[m++]; - nz[j][2] += buf[m++]; - - } -} - -/* ---------------------------------------------------------------------- */ - void *PairOxdnaHbond::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"nx") == 0) return (void *) nx; - if (strcmp(str,"ny") == 0) return (void *) ny; - if (strcmp(str,"nz") == 0) return (void *) nz; - if (strcmp(str,"epsilon_hb") == 0) return (void *) epsilon_hb; if (strcmp(str,"a_hb") == 0) return (void *) a_hb; if (strcmp(str,"cut_hb_0") == 0) return (void *) cut_hb_0; diff --git a/src/CG-DNA/pair_oxdna_stk.cpp b/src/CG-DNA/pair_oxdna_stk.cpp index ab9152a8d7..f156062463 100644 --- a/src/CG-DNA/pair_oxdna_stk.cpp +++ b/src/CG-DNA/pair_oxdna_stk.cpp @@ -229,10 +229,9 @@ void PairOxdnaStk::compute(int eflag, int vflag) // vectors COM-backbone site, COM-stacking site in lab frame double ra_cs[3],ra_cst[3]; double rb_cs[3],rb_cst[3]; - - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; + // Cartesian unit vectors in lab frame + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -258,6 +257,12 @@ void PairOxdnaStk::compute(int eflag, int vflag) evdwl = 0.0; ev_init(eflag,vflag); + + // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + ny_xtrct = (double **) force->pair->extract("ny",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); // loop over stacking interaction neighbors using bond topology @@ -277,10 +282,13 @@ void PairOxdnaStk::compute(int eflag, int vflag) // 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); + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + //(a/b)y/z not needed here as oxDNA(1) co-linear // vector COM a - stacking site a ra_cst[0] = d_cst*ax[0]; @@ -337,6 +345,13 @@ void PairOxdnaStk::compute(int eflag, int vflag) // early rejection criterium if (f1) { + + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; // theta4 angle and correction cost4 = MathExtra::dot3(bz,az); @@ -361,6 +376,13 @@ void PairOxdnaStk::compute(int eflag, int vflag) // early rejection criterium if (f4t5) { + + ay[0] = ny_xtrct[a][0]; + ay[1] = ny_xtrct[a][1]; + ay[2] = ny_xtrct[a][2]; + by[0] = ny_xtrct[b][0]; + by[1] = ny_xtrct[b][1]; + by[2] = ny_xtrct[b][2]; cost6p = MathExtra::dot3(delr_st_norm,az); if (cost6p > 1.0) cost6p = 1.0; diff --git a/src/CG-DNA/pair_oxdna_stk.h b/src/CG-DNA/pair_oxdna_stk.h index 8695f8fd36..ada5cdcc6a 100644 --- a/src/CG-DNA/pair_oxdna_stk.h +++ b/src/CG-DNA/pair_oxdna_stk.h @@ -59,6 +59,9 @@ class PairOxdnaStk : public Pair { double **b_st6, **dtheta_st6_c; double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c; double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c; + + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage + int seqdepflag; diff --git a/src/CG-DNA/pair_oxdna_xstk.cpp b/src/CG-DNA/pair_oxdna_xstk.cpp index 84e67666ee..046659e5c6 100644 --- a/src/CG-DNA/pair_oxdna_xstk.cpp +++ b/src/CG-DNA/pair_oxdna_xstk.cpp @@ -128,10 +128,9 @@ void PairOxdnaXstk::compute(int eflag, int vflag) double d_chb=+0.4; // vectors COM-h-bonding site in lab frame double ra_chb[3],rb_chb[3]; - - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; + // Cartesian unit vectors in lab frame + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -159,6 +158,12 @@ void PairOxdnaXstk::compute(int eflag, int vflag) alist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; + + // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + ny_xtrct = (double **) force->pair->extract("ny",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); // loop over pair interaction neighbors of my atoms @@ -167,8 +172,10 @@ void PairOxdnaXstk::compute(int eflag, int vflag) a = alist[ia]; atype = type[a]; - qa=bonus[ellipsoid[a]].quat; - MathExtra::q_to_exyz(qa,ax,ay,az); + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + //a(y/z) not needed here as oxDNA(1) co-linear ra_chb[0] = d_chb*ax[0]; ra_chb[1] = d_chb*ax[1]; @@ -185,8 +192,10 @@ void PairOxdnaXstk::compute(int eflag, int vflag) btype = type[b]; - qb=bonus[ellipsoid[b]].quat; - MathExtra::q_to_exyz(qb,bx,by,bz); + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + //b(y/z) not needed here as oxDNA(1) co-linear rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; @@ -244,6 +253,13 @@ void PairOxdnaXstk::compute(int eflag, int vflag) // early rejection criterium if (f4t3) { + + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; diff --git a/src/CG-DNA/pair_oxdna_xstk.h b/src/CG-DNA/pair_oxdna_xstk.h index 30089f53fb..c450755356 100644 --- a/src/CG-DNA/pair_oxdna_xstk.h +++ b/src/CG-DNA/pair_oxdna_xstk.h @@ -65,6 +65,9 @@ class PairOxdnaXstk : public Pair { double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast; double **b_xst8, **dtheta_xst8_c; + + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage + virtual void allocate(); }; diff --git a/src/atom.cpp b/src/atom.cpp index d368ce7f85..8bf361ab95 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -650,7 +650,6 @@ void Atom::set_atomflag_defaults() rho_flag = esph_flag = cv_flag = vest_flag = 0; dpd_flag = edpd_flag = tdpd_flag = 0; sp_flag = 0; - lrefpos_flag = 0; x0_flag = 0; smd_flag = damage_flag = 0; mesont_flag = 0; diff --git a/src/atom.h b/src/atom.h index b37e0e49e0..25a56de1c4 100644 --- a/src/atom.h +++ b/src/atom.h @@ -190,10 +190,6 @@ class Atom : protected Pointers { // SPIN package int sp_flag; - - // CG-DNA package - - int lrefpos_flag; // MACHDYN package From b8970366e01146934e99f1ae9e60d41000c1d203 Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Wed, 3 Nov 2021 10:59:04 +0000 Subject: [PATCH 04/13] Fixed whitespace issues --- src/CG-DNA/bond_oxdna_fene.cpp | 4 +- src/CG-DNA/bond_oxdna_fene.h | 2 +- src/CG-DNA/pair_oxdna_coaxstk.cpp | 37 +++++++++---------- src/CG-DNA/pair_oxdna_coaxstk.h | 2 +- src/CG-DNA/pair_oxdna_excv.cpp | 61 +++++++++++++++---------------- src/CG-DNA/pair_oxdna_excv.h | 3 +- src/CG-DNA/pair_oxdna_hbond.cpp | 34 ++++++++--------- src/CG-DNA/pair_oxdna_hbond.h | 2 +- src/CG-DNA/pair_oxdna_stk.cpp | 39 +++++++++----------- src/CG-DNA/pair_oxdna_stk.h | 2 +- src/CG-DNA/pair_oxdna_xstk.cpp | 27 ++++++-------- src/CG-DNA/pair_oxdna_xstk.h | 2 +- 12 files changed, 103 insertions(+), 112 deletions(-) diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp index c46556768c..317b6b3c9b 100644 --- a/src/CG-DNA/bond_oxdna_fene.cpp +++ b/src/CG-DNA/bond_oxdna_fene.cpp @@ -174,7 +174,7 @@ void BondOxdnaFene::compute(int eflag, int vflag) ebond = 0.0; ev_init(eflag, vflag); - // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -194,7 +194,7 @@ void BondOxdnaFene::compute(int eflag, int vflag) bx[0] = nx_xtrct[b][0]; bx[1] = nx_xtrct[b][1]; bx[2] = nx_xtrct[b][2]; - //(a/b)y/z not needed here as oxDNA(1) co-linear + // (a/b)y/z not needed here as oxDNA(1) co-linear // vector COM-backbone site a and b compute_interaction_sites(ax, ay, az, ra_cs); diff --git a/src/CG-DNA/bond_oxdna_fene.h b/src/CG-DNA/bond_oxdna_fene.h index b24757781f..e6020892e4 100644 --- a/src/CG-DNA/bond_oxdna_fene.h +++ b/src/CG-DNA/bond_oxdna_fene.h @@ -40,7 +40,7 @@ class BondOxdnaFene : public Bond { protected: double *k, *Delta, *r0; // FENE - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors void allocate(); void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double); diff --git a/src/CG-DNA/pair_oxdna_coaxstk.cpp b/src/CG-DNA/pair_oxdna_coaxstk.cpp index f76ad33147..a5b94403da 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna_coaxstk.cpp @@ -106,9 +106,6 @@ PairOxdnaCoaxstk::~PairOxdnaCoaxstk() void PairOxdnaCoaxstk::compute(int eflag, int vflag) { - - //printf("\n Coaxstk HERE, proc = %d \n", comm->me); - double delf[3],delt[3],delta[3],deltb[3]; // force, torque increment; double evdwl,fpair,finc,tpair,factor_lj; double v1tmp[3],v2tmp[3],v3tmp[3]; @@ -162,7 +159,7 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; - // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -176,9 +173,9 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) atype = type[a]; ax[0] = nx_xtrct[a][0]; - ax[1] = nx_xtrct[a][1]; - ax[2] = nx_xtrct[a][2]; - //a(y/z) not needed here as oxDNA(1) co-linear + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + // a(y/z) not needed here as oxDNA(1) co-linear // vector COM a - stacking site a ra_cst[0] = d_cst*ax[0]; @@ -202,9 +199,9 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) btype = type[b]; bx[0] = nx_xtrct[b][0]; - bx[1] = nx_xtrct[b][1]; - bx[2] = nx_xtrct[b][2]; - //b(y/z) not needed here as oxDNA(1) co-linear + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + // b(y/z) not needed here as oxDNA(1) co-linear // vector COM b - stacking site b rb_cst[0] = d_cst*bx[0]; @@ -256,12 +253,12 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) // early rejection criterium if (f4t1) { - az[0] = nz_xtrct[a][0]; - az[1] = nz_xtrct[a][1]; - az[2] = nz_xtrct[a][2]; - bz[0] = nz_xtrct[b][0]; - bz[1] = nz_xtrct[b][1]; - bz[2] = nz_xtrct[b][2]; + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; @@ -398,9 +395,9 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) // cosphi3 and cosphi4 (=cosphi3) force and virial if (cosphi3) { - ay[0] = ny_xtrct[a][0]; - ay[1] = ny_xtrct[a][1]; - ay[2] = ny_xtrct[a][2]; + ay[0] = ny_xtrct[a][0]; + ay[1] = ny_xtrct[a][1]; + ay[2] = ny_xtrct[a][2]; finc = -f2 * f4t1* f4t4 * f4t5 * f4t6 * 2.0 * f5c3 * df5c3 * factor_lj; fpair += finc; @@ -408,6 +405,7 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) gamma = d_cs - d_cst; gammacub = gamma * gamma * gamma; rinv_ss_cub = rinv_ss * rinv_ss * rinv_ss; + aybx = MathExtra::dot3(ay,bx); azbx = MathExtra::dot3(az,bx); rax = MathExtra::dot3(delr_st_norm,ax); @@ -547,6 +545,7 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) gamma = d_cs - d_cst; gammacub = gamma * gamma * gamma; rinv_ss_cub = rinv_ss * rinv_ss * rinv_ss; + aybx = MathExtra::dot3(ay,bx); azbx = MathExtra::dot3(az,bx); rax = MathExtra::dot3(delr_st_norm,ax); diff --git a/src/CG-DNA/pair_oxdna_coaxstk.h b/src/CG-DNA/pair_oxdna_coaxstk.h index d634ccdf00..050fd5b06a 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.h +++ b/src/CG-DNA/pair_oxdna_coaxstk.h @@ -63,7 +63,7 @@ class PairOxdnaCoaxstk : public Pair { double **a_cxst3p, **cosphi_cxst3p_ast, **b_cxst3p, **cosphi_cxst3p_c; double **a_cxst4p, **cosphi_cxst4p_ast, **b_cxst4p, **cosphi_cxst4p_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxdna_excv.cpp b/src/CG-DNA/pair_oxdna_excv.cpp index 2edad16c83..096d9b5171 100644 --- a/src/CG-DNA/pair_oxdna_excv.cpp +++ b/src/CG-DNA/pair_oxdna_excv.cpp @@ -50,7 +50,7 @@ PairOxdnaExcv::~PairOxdnaExcv() { if (allocated) { - memory->destroy(nx); + memory->destroy(nx); memory->destroy(ny); memory->destroy(nz); @@ -156,24 +156,23 @@ void PairOxdnaExcv::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; - // loop over all local atoms, handle calculation of local reference frame - + // loop over all local atoms, calculation of local reference frame for (in = 0; in < atom->nlocal; in++) { int n = alist[in]; - double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame - qn=bonus[ellipsoid[n]].quat; - MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp); - - nx[n][0] = nx_temp[0]; - nx[n][1] = nx_temp[1]; - nx[n][2] = nx_temp[2]; - ny[n][0] = ny_temp[0]; - ny[n][1] = ny_temp[1]; - ny[n][2] = ny_temp[2]; - nz[n][0] = nz_temp[0]; - nz[n][1] = nz_temp[1]; + + qn=bonus[ellipsoid[n]].quat; + MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp); + + nx[n][0] = nx_temp[0]; + nx[n][1] = nx_temp[1]; + nx[n][2] = nx_temp[2]; + ny[n][0] = ny_temp[0]; + ny[n][1] = ny_temp[1]; + ny[n][2] = ny_temp[2]; + nz[n][0] = nz_temp[0]; + nz[n][1] = nz_temp[1]; nz[n][2] = nz_temp[2]; } @@ -188,14 +187,14 @@ void PairOxdnaExcv::compute(int eflag, int vflag) atype = type[a]; ax[0] = nx[a][0]; - ax[1] = nx[a][1]; - ax[2] = nx[a][2]; - ay[0] = ny[a][0]; - ay[1] = ny[a][1]; - ay[2] = ny[a][2]; - az[0] = nz[a][0]; - az[1] = nz[a][1]; - az[2] = nz[a][2]; + ax[1] = nx[a][1]; + ax[2] = nx[a][2]; + ay[0] = ny[a][0]; + ay[1] = ny[a][1]; + ay[2] = ny[a][2]; + az[0] = nz[a][0]; + az[1] = nz[a][1]; + az[2] = nz[a][2]; // vector COM - backbone and base site a compute_interaction_sites(ax,ay,az,ra_cs,ra_cb); @@ -220,14 +219,14 @@ void PairOxdnaExcv::compute(int eflag, int vflag) btype = type[b]; bx[0] = nx[b][0]; - bx[1] = nx[b][1]; - bx[2] = nx[b][2]; - by[0] = ny[b][0]; - by[1] = ny[b][1]; - by[2] = ny[b][2]; - bz[0] = nz[b][0]; - bz[1] = nz[b][1]; - bz[2] = nz[b][2]; + bx[1] = nx[b][1]; + bx[2] = nx[b][2]; + by[0] = ny[b][0]; + by[1] = ny[b][1]; + by[2] = ny[b][2]; + bz[0] = nz[b][0]; + bz[1] = nz[b][1]; + bz[2] = nz[b][2]; // vector COM - backbone and base site b compute_interaction_sites(bx,by,bz,rb_cs,rb_cb); diff --git a/src/CG-DNA/pair_oxdna_excv.h b/src/CG-DNA/pair_oxdna_excv.h index 1d2a871dd9..9abd97fb3d 100644 --- a/src/CG-DNA/pair_oxdna_excv.h +++ b/src/CG-DNA/pair_oxdna_excv.h @@ -56,8 +56,7 @@ class PairOxdnaExcv : public Pair { double **epsilon_bb, **sigma_bb, **cut_bb_ast, **cutsq_bb_ast; double **lj1_bb, **lj2_bb, **b_bb, **cut_bb_c, **cutsq_bb_c; - // per-atom arrays for q_to_exyz storage - double **nx, **ny, **nz; + double **nx, **ny, **nz; // per-atom arrays for local unit vectors virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxdna_hbond.cpp b/src/CG-DNA/pair_oxdna_hbond.cpp index f4baa74787..07a1ce332b 100644 --- a/src/CG-DNA/pair_oxdna_hbond.cpp +++ b/src/CG-DNA/pair_oxdna_hbond.cpp @@ -179,7 +179,7 @@ void PairOxdnaHbond::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; - // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -193,10 +193,10 @@ void PairOxdnaHbond::compute(int eflag, int vflag) atype = type[a]; ax[0] = nx_xtrct[a][0]; - ax[1] = nx_xtrct[a][1]; - ax[2] = nx_xtrct[a][2]; - - ra_chb[0] = d_chb*ax[0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + + ra_chb[0] = d_chb*ax[0]; ra_chb[1] = d_chb*ax[1]; ra_chb[2] = d_chb*ax[2]; @@ -211,11 +211,11 @@ void PairOxdnaHbond::compute(int eflag, int vflag) btype = type[b]; - bx[0] = nx_xtrct[b][0]; - bx[1] = nx_xtrct[b][1]; - bx[2] = nx_xtrct[b][2]; - - rb_chb[0] = d_chb*bx[0]; + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + + rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; rb_chb[2] = d_chb*bx[2]; @@ -272,12 +272,12 @@ void PairOxdnaHbond::compute(int eflag, int vflag) // early rejection criterium if (f4t3) { - az[0] = nz_xtrct[a][0]; - az[1] = nz_xtrct[a][1]; - az[2] = nz_xtrct[a][2]; - bz[0] = nz_xtrct[b][0]; - bz[1] = nz_xtrct[b][1]; - bz[2] = nz_xtrct[b][2]; + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; @@ -1249,4 +1249,4 @@ void *PairOxdnaHbond::extract(const char *str, int &dim) if (strcmp(str,"dtheta_hb8_c") == 0) return (void *) dtheta_hb8_c; return nullptr; -} \ No newline at end of file +} diff --git a/src/CG-DNA/pair_oxdna_hbond.h b/src/CG-DNA/pair_oxdna_hbond.h index 90ede130ba..1a138e82de 100644 --- a/src/CG-DNA/pair_oxdna_hbond.h +++ b/src/CG-DNA/pair_oxdna_hbond.h @@ -67,7 +67,7 @@ class PairOxdnaHbond : public Pair { double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast; double **b_hb8, **dtheta_hb8_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors int seqdepflag; diff --git a/src/CG-DNA/pair_oxdna_stk.cpp b/src/CG-DNA/pair_oxdna_stk.cpp index f156062463..7665abbb79 100644 --- a/src/CG-DNA/pair_oxdna_stk.cpp +++ b/src/CG-DNA/pair_oxdna_stk.cpp @@ -212,9 +212,6 @@ void PairOxdnaStk::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, void PairOxdnaStk::compute(int eflag, int vflag) { - - //printf("\n Stk HERE, proc = %d \n", comm->me); - double delf[3],delta[3],deltb[3]; // force, torque increment; double evdwl,fpair,finc,tpair; double delr_ss[3],delr_ss_norm[3],rsq_ss,r_ss,rinv_ss; @@ -258,7 +255,7 @@ void PairOxdnaStk::compute(int eflag, int vflag) evdwl = 0.0; ev_init(eflag,vflag); - // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -282,13 +279,13 @@ void PairOxdnaStk::compute(int eflag, int vflag) // a now in 3' direction, b in 5' direction - ax[0] = nx_xtrct[a][0]; - ax[1] = nx_xtrct[a][1]; - ax[2] = nx_xtrct[a][2]; - bx[0] = nx_xtrct[b][0]; - bx[1] = nx_xtrct[b][1]; - bx[2] = nx_xtrct[b][2]; - //(a/b)y/z not needed here as oxDNA(1) co-linear + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + // (a/b)y/z not needed here as oxDNA(1) co-linear // vector COM a - stacking site a ra_cst[0] = d_cst*ax[0]; @@ -346,12 +343,12 @@ void PairOxdnaStk::compute(int eflag, int vflag) // early rejection criterium if (f1) { - az[0] = nz_xtrct[a][0]; - az[1] = nz_xtrct[a][1]; - az[2] = nz_xtrct[a][2]; - bz[0] = nz_xtrct[b][0]; - bz[1] = nz_xtrct[b][1]; - bz[2] = nz_xtrct[b][2]; + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; // theta4 angle and correction cost4 = MathExtra::dot3(bz,az); @@ -377,12 +374,12 @@ void PairOxdnaStk::compute(int eflag, int vflag) // early rejection criterium if (f4t5) { - ay[0] = ny_xtrct[a][0]; + ay[0] = ny_xtrct[a][0]; ay[1] = ny_xtrct[a][1]; - ay[2] = ny_xtrct[a][2]; - by[0] = ny_xtrct[b][0]; + ay[2] = ny_xtrct[a][2]; + by[0] = ny_xtrct[b][0]; by[1] = ny_xtrct[b][1]; - by[2] = ny_xtrct[b][2]; + by[2] = ny_xtrct[b][2]; cost6p = MathExtra::dot3(delr_st_norm,az); if (cost6p > 1.0) cost6p = 1.0; diff --git a/src/CG-DNA/pair_oxdna_stk.h b/src/CG-DNA/pair_oxdna_stk.h index ada5cdcc6a..1f8a3a9249 100644 --- a/src/CG-DNA/pair_oxdna_stk.h +++ b/src/CG-DNA/pair_oxdna_stk.h @@ -60,7 +60,7 @@ class PairOxdnaStk : public Pair { double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c; double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors int seqdepflag; diff --git a/src/CG-DNA/pair_oxdna_xstk.cpp b/src/CG-DNA/pair_oxdna_xstk.cpp index 046659e5c6..bfb0d185a9 100644 --- a/src/CG-DNA/pair_oxdna_xstk.cpp +++ b/src/CG-DNA/pair_oxdna_xstk.cpp @@ -111,9 +111,6 @@ PairOxdnaXstk::~PairOxdnaXstk() void PairOxdnaXstk::compute(int eflag, int vflag) { - - //printf("\n Xstk HERE, proc = %d \n", comm->me); - double delf[3],delta[3],deltb[3]; // force, torque increment; double evdwl,fpair,finc,tpair,factor_lj; double delr_hb[3],delr_hb_norm[3],rsq_hb,r_hb,rinv_hb; @@ -159,7 +156,7 @@ void PairOxdnaXstk::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; - // n(x/y/z)_xtrct = extracted q_to_exyz from oxdna_excv + // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -173,9 +170,9 @@ void PairOxdnaXstk::compute(int eflag, int vflag) atype = type[a]; ax[0] = nx_xtrct[a][0]; - ax[1] = nx_xtrct[a][1]; - ax[2] = nx_xtrct[a][2]; - //a(y/z) not needed here as oxDNA(1) co-linear + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + // a(y/z) not needed here as oxDNA(1) co-linear ra_chb[0] = d_chb*ax[0]; ra_chb[1] = d_chb*ax[1]; @@ -193,9 +190,9 @@ void PairOxdnaXstk::compute(int eflag, int vflag) btype = type[b]; bx[0] = nx_xtrct[b][0]; - bx[1] = nx_xtrct[b][1]; - bx[2] = nx_xtrct[b][2]; - //b(y/z) not needed here as oxDNA(1) co-linear + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + // b(y/z) not needed here as oxDNA(1) co-linear rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; @@ -254,12 +251,12 @@ void PairOxdnaXstk::compute(int eflag, int vflag) // early rejection criterium if (f4t3) { - az[0] = nz_xtrct[a][0]; + az[0] = nz_xtrct[a][0]; az[1] = nz_xtrct[a][1]; - az[2] = nz_xtrct[a][2]; - bz[0] = nz_xtrct[b][0]; - bz[1] = nz_xtrct[b][1]; - bz[2] = nz_xtrct[b][2]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; diff --git a/src/CG-DNA/pair_oxdna_xstk.h b/src/CG-DNA/pair_oxdna_xstk.h index c450755356..2ad4c6027c 100644 --- a/src/CG-DNA/pair_oxdna_xstk.h +++ b/src/CG-DNA/pair_oxdna_xstk.h @@ -66,7 +66,7 @@ class PairOxdnaXstk : public Pair { double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast; double **b_xst8, **dtheta_xst8_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for q_to_exyz storage + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors virtual void allocate(); From fc4fdd09efcd7bda420bae5fc4eeefbbae4134c8 Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Wed, 3 Nov 2021 11:44:29 +0000 Subject: [PATCH 05/13] Fixed more whitespace issues --- src/CG-DNA/bond_oxdna_fene.cpp | 29 ++++++++++++++--------------- src/CG-DNA/pair_oxdna_coaxstk.cpp | 2 -- src/CG-DNA/pair_oxdna_excv.cpp | 3 --- src/CG-DNA/pair_oxdna_hbond.cpp | 2 +- 4 files changed, 15 insertions(+), 21 deletions(-) diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp index 317b6b3c9b..a8cd52e651 100644 --- a/src/CG-DNA/bond_oxdna_fene.cpp +++ b/src/CG-DNA/bond_oxdna_fene.cpp @@ -146,14 +146,13 @@ void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, doub ------------------------------------------------------------------------- */ void BondOxdnaFene::compute(int eflag, int vflag) { - - int a, b, in, type; - double delf[3], delta[3], deltb[3]; // force, torque increment;; - double delr[3], ebond, fbond; - double rsq, Deltasq, rlogarg; - double r, rr0, rr0sq; + int a,b,in,type; + double delf[3],delta[3],deltb[3]; // force, torque increment;; + double delr[3],ebond,fbond; + double rsq,Deltasq,rlogarg; + double r,rr0,rr0sq; // vectors COM-backbone site in lab frame - double ra_cs[3], rb_cs[3]; + double ra_cs[3],rb_cs[3]; // Cartesian unit vectors in lab frame double ax[3],ay[3],az[3]; double bx[3],by[3],bz[3]; @@ -188,13 +187,13 @@ void BondOxdnaFene::compute(int eflag, int vflag) b = bondlist[in][0]; type = bondlist[in][2]; - ax[0] = nx_xtrct[a][0]; - ax[1] = nx_xtrct[a][1]; - ax[2] = nx_xtrct[a][2]; - bx[0] = nx_xtrct[b][0]; - bx[1] = nx_xtrct[b][1]; - bx[2] = nx_xtrct[b][2]; - // (a/b)y/z not needed here as oxDNA(1) co-linear + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + // (a/b)y/z not needed here as oxDNA(1) co-linear // vector COM-backbone site a and b compute_interaction_sites(ax, ay, az, ra_cs); @@ -266,7 +265,7 @@ void BondOxdnaFene::compute(int eflag, int vflag) if (evflag) ev_tally_xyz(a, b, nlocal, newton_bond, ebond, delf[0], delf[1], delf[2], x[a][0] - x[b][0], x[a][1] - x[b][1], x[a][2] - x[b][2]); - } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/CG-DNA/pair_oxdna_coaxstk.cpp b/src/CG-DNA/pair_oxdna_coaxstk.cpp index a5b94403da..e03de69f57 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna_coaxstk.cpp @@ -405,7 +405,6 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) gamma = d_cs - d_cst; gammacub = gamma * gamma * gamma; rinv_ss_cub = rinv_ss * rinv_ss * rinv_ss; - aybx = MathExtra::dot3(ay,bx); azbx = MathExtra::dot3(az,bx); rax = MathExtra::dot3(delr_st_norm,ax); @@ -545,7 +544,6 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) gamma = d_cs - d_cst; gammacub = gamma * gamma * gamma; rinv_ss_cub = rinv_ss * rinv_ss * rinv_ss; - aybx = MathExtra::dot3(ay,bx); azbx = MathExtra::dot3(az,bx); rax = MathExtra::dot3(delr_st_norm,ax); diff --git a/src/CG-DNA/pair_oxdna_excv.cpp b/src/CG-DNA/pair_oxdna_excv.cpp index 096d9b5171..2ae5995ebb 100644 --- a/src/CG-DNA/pair_oxdna_excv.cpp +++ b/src/CG-DNA/pair_oxdna_excv.cpp @@ -115,9 +115,6 @@ void PairOxdnaExcv::compute_interaction_sites(double e1[3], double /*e2*/[3], void PairOxdnaExcv::compute(int eflag, int vflag) { - - //printf("\n ExcVol HERE, proc = %d \n", comm->me); - double delf[3],delta[3],deltb[3]; // force, torque increment; double evdwl,fpair,factor_lj; double rtmp_s[3],rtmp_b[3]; diff --git a/src/CG-DNA/pair_oxdna_hbond.cpp b/src/CG-DNA/pair_oxdna_hbond.cpp index 07a1ce332b..eb04d74f2e 100644 --- a/src/CG-DNA/pair_oxdna_hbond.cpp +++ b/src/CG-DNA/pair_oxdna_hbond.cpp @@ -218,7 +218,7 @@ void PairOxdnaHbond::compute(int eflag, int vflag) rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; rb_chb[2] = d_chb*bx[2]; - + // vector h-bonding site b to a delr_hb[0] = x[a][0] + ra_chb[0] - x[b][0] - rb_chb[0]; delr_hb[1] = x[a][1] + ra_chb[1] - x[b][1] - rb_chb[1]; From e8f6024eae460ca3fb6dbff7508f0b873293304f Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Wed, 3 Nov 2021 11:46:56 +0000 Subject: [PATCH 06/13] Fixed more whitespace issues --- src/CG-DNA/bond_oxdna_fene.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp index a8cd52e651..3652856a52 100644 --- a/src/CG-DNA/bond_oxdna_fene.cpp +++ b/src/CG-DNA/bond_oxdna_fene.cpp @@ -145,7 +145,7 @@ void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, doub s=sugar-phosphate backbone site, b=base site, st=stacking site ------------------------------------------------------------------------- */ void BondOxdnaFene::compute(int eflag, int vflag) -{ +{ int a,b,in,type; double delf[3],delta[3],deltb[3]; // force, torque increment;; double delr[3],ebond,fbond; From f733453f0547ecb206e5077f168c313a7a55194a Mon Sep 17 00:00:00 2001 From: Lewis M Russell <71833891+lrussell676@users.noreply.github.com> Date: Wed, 10 Nov 2021 09:25:13 +0000 Subject: [PATCH 07/13] Intpos (#11) * hbond comm added for rsq_hb * lrefpos removed, extract scaled for oxDNA1 * Update pair_oxdna_hbond.cpp * nxyz extract scaled across DNA2/RNA2 * oxDNA2/RNA2 updated to match oxDNA styling from upstream merge * whitespace corrections also removed unnecessary local unit vector from oxRNA2_xstk * whitespace correction in oxdna_coaxstk --- src/CG-DNA/bond_oxdna_fene.cpp | 17 +++++++++++--- src/CG-DNA/pair_oxdna2_coaxstk.cpp | 30 ++++++++++++++++++------- src/CG-DNA/pair_oxdna2_coaxstk.h | 3 +++ src/CG-DNA/pair_oxdna2_dh.cpp | 36 +++++++++++++++++++++++------- src/CG-DNA/pair_oxdna2_dh.h | 3 +++ src/CG-DNA/pair_oxdna_coaxstk.cpp | 10 ++++----- src/CG-DNA/pair_oxdna_coaxstk.h | 2 +- src/CG-DNA/pair_oxdna_excv.cpp | 26 ++++++++++----------- src/CG-DNA/pair_oxdna_hbond.cpp | 12 +++++----- src/CG-DNA/pair_oxdna_hbond.h | 4 ++-- src/CG-DNA/pair_oxdna_stk.cpp | 12 +++++----- src/CG-DNA/pair_oxdna_stk.h | 2 +- src/CG-DNA/pair_oxdna_xstk.cpp | 6 ++--- src/CG-DNA/pair_oxdna_xstk.h | 2 +- src/CG-DNA/pair_oxrna2_stk.cpp | 36 +++++++++++++++++++++++------- src/CG-DNA/pair_oxrna2_stk.h | 3 +++ src/CG-DNA/pair_oxrna2_xstk.cpp | 29 ++++++++++++++++++------ src/CG-DNA/pair_oxrna2_xstk.h | 3 +++ 18 files changed, 164 insertions(+), 72 deletions(-) diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp index 3652856a52..4db0fb4cf1 100644 --- a/src/CG-DNA/bond_oxdna_fene.cpp +++ b/src/CG-DNA/bond_oxdna_fene.cpp @@ -172,8 +172,8 @@ void BondOxdnaFene::compute(int eflag, int vflag) ebond = 0.0; ev_init(eflag, vflag); - - // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv + + // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -190,10 +190,21 @@ void BondOxdnaFene::compute(int eflag, int vflag) ax[0] = nx_xtrct[a][0]; ax[1] = nx_xtrct[a][1]; ax[2] = nx_xtrct[a][2]; + ay[0] = ny_xtrct[a][0]; + ay[1] = ny_xtrct[a][1]; + ay[2] = ny_xtrct[a][2]; + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; bx[0] = nx_xtrct[b][0]; bx[1] = nx_xtrct[b][1]; bx[2] = nx_xtrct[b][2]; - // (a/b)y/z not needed here as oxDNA(1) co-linear + by[0] = ny_xtrct[b][0]; + by[1] = ny_xtrct[b][1]; + by[2] = ny_xtrct[b][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; // vector COM-backbone site a and b compute_interaction_sites(ax, ay, az, ra_cs); diff --git a/src/CG-DNA/pair_oxdna2_coaxstk.cpp b/src/CG-DNA/pair_oxdna2_coaxstk.cpp index e8ab197be1..8b3609983b 100644 --- a/src/CG-DNA/pair_oxdna2_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna2_coaxstk.cpp @@ -118,9 +118,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag) double ra_cs[3],ra_cst[3]; double rb_cs[3],rb_cst[3]; - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; + // Cartesian unit vectors in lab frame + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -149,6 +149,11 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // n(x/z)_xtrct = extracted local unit vectors from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); + // loop over pair interaction neighbors of my atoms for (ia = 0; ia < anum; ia++) { @@ -156,8 +161,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag) a = alist[ia]; atype = type[a]; - qa=bonus[ellipsoid[a]].quat; - MathExtra::q_to_exyz(qa,ax,ay,az); + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; // vector COM a - stacking site a ra_cst[0] = d_cst*ax[0]; @@ -180,8 +186,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag) btype = type[b]; - qb=bonus[ellipsoid[b]].quat; - MathExtra::q_to_exyz(qb,bx,by,bz); + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; // vector COM b - stacking site b rb_cst[0] = d_cst*bx[0]; @@ -232,6 +239,13 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag) // early rejection criterium if (f4f6t1) { + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; + cost4 = MathExtra::dot3(az,bz); if (cost4 > 1.0) cost4 = 1.0; if (cost4 < -1.0) cost4 = -1.0; @@ -306,7 +320,7 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag) DF4(theta6p, a_cxst6[atype][btype], theta_cxst6_0[atype][btype], dtheta_cxst6_ast[atype][btype], b_cxst6[atype][btype], dtheta_cxst6_c[atype][btype])*rsint; - // force, torque and virial contribution for forces between stacking sites + // force, torque and virial contribution for forces between stacking sites fpair = 0.0; diff --git a/src/CG-DNA/pair_oxdna2_coaxstk.h b/src/CG-DNA/pair_oxdna2_coaxstk.h index 2a323be579..7b376fa368 100644 --- a/src/CG-DNA/pair_oxdna2_coaxstk.h +++ b/src/CG-DNA/pair_oxdna2_coaxstk.h @@ -61,6 +61,9 @@ class PairOxdna2Coaxstk : public Pair { double **AA_cxst1, **BB_cxst1; + // per-atom arrays for local unit vectors + double **nx_xtrct, **nz_xtrct; + virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxdna2_dh.cpp b/src/CG-DNA/pair_oxdna2_dh.cpp index 1a587e82e1..cf79fd07cf 100644 --- a/src/CG-DNA/pair_oxdna2_dh.cpp +++ b/src/CG-DNA/pair_oxdna2_dh.cpp @@ -87,10 +87,9 @@ void PairOxdna2Dh::compute(int eflag, int vflag) // vectors COM-backbone sites in lab frame double ra_cs[3],rb_cs[3]; - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; - double *special_lj = force->special_lj; + // Cartesian unit vectors in lab frame + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -100,6 +99,7 @@ void PairOxdna2Dh::compute(int eflag, int vflag) int nlocal = atom->nlocal; int newton_pair = force->newton_pair; int *alist,*blist,*numneigh,**firstneigh; + double *special_lj = force->special_lj; AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); AtomVecEllipsoid::Bonus *bonus = avec->bonus; @@ -115,6 +115,12 @@ void PairOxdna2Dh::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + ny_xtrct = (double **) force->pair->extract("ny",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); + // loop over pair interaction neighbors of my atoms for (ia = 0; ia < anum; ia++) { @@ -122,8 +128,15 @@ void PairOxdna2Dh::compute(int eflag, int vflag) a = alist[ia]; atype = type[a]; - qa=bonus[ellipsoid[a]].quat; - MathExtra::q_to_exyz(qa,ax,ay,az); + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + ay[0] = ny_xtrct[a][0]; + ay[1] = ny_xtrct[a][1]; + ay[2] = ny_xtrct[a][2]; + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; // vector COM-backbone site a compute_interaction_sites(ax,ay,az,ra_cs); @@ -142,8 +155,15 @@ void PairOxdna2Dh::compute(int eflag, int vflag) b &= NEIGHMASK; btype = type[b]; - qb=bonus[ellipsoid[b]].quat; - MathExtra::q_to_exyz(qb,bx,by,bz); + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + by[0] = ny_xtrct[b][0]; + by[1] = ny_xtrct[b][1]; + by[2] = ny_xtrct[b][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; // vector COM-backbone site b compute_interaction_sites(bx,by,bz,rb_cs); diff --git a/src/CG-DNA/pair_oxdna2_dh.h b/src/CG-DNA/pair_oxdna2_dh.h index 925b53e322..d7075f2988 100644 --- a/src/CG-DNA/pair_oxdna2_dh.h +++ b/src/CG-DNA/pair_oxdna2_dh.h @@ -46,6 +46,9 @@ class PairOxdna2Dh : public Pair { double **qeff_dh_pf, **kappa_dh; double **b_dh, **cut_dh_ast, **cutsq_dh_ast, **cut_dh_c, **cutsq_dh_c; + // per-atom arrays for local unit vectors + double **nx_xtrct, **ny_xtrct, **nz_xtrct; + virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxdna_coaxstk.cpp b/src/CG-DNA/pair_oxdna_coaxstk.cpp index e03de69f57..e6cf59f072 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna_coaxstk.cpp @@ -159,7 +159,7 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; - // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv + // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -252,7 +252,7 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) // early rejection criterium if (f4t1) { - + az[0] = nz_xtrct[a][0]; az[1] = nz_xtrct[a][1]; az[2] = nz_xtrct[a][2]; @@ -395,9 +395,9 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) // cosphi3 and cosphi4 (=cosphi3) force and virial if (cosphi3) { - ay[0] = ny_xtrct[a][0]; - ay[1] = ny_xtrct[a][1]; - ay[2] = ny_xtrct[a][2]; + ay[0] = ny_xtrct[a][0]; + ay[1] = ny_xtrct[a][1]; + ay[2] = ny_xtrct[a][2]; finc = -f2 * f4t1* f4t4 * f4t5 * f4t6 * 2.0 * f5c3 * df5c3 * factor_lj; fpair += finc; diff --git a/src/CG-DNA/pair_oxdna_coaxstk.h b/src/CG-DNA/pair_oxdna_coaxstk.h index 050fd5b06a..28813a034f 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.h +++ b/src/CG-DNA/pair_oxdna_coaxstk.h @@ -62,7 +62,7 @@ class PairOxdnaCoaxstk : public Pair { double **a_cxst3p, **cosphi_cxst3p_ast, **b_cxst3p, **cosphi_cxst3p_c; double **a_cxst4p, **cosphi_cxst4p_ast, **b_cxst4p, **cosphi_cxst4p_c; - + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors virtual void allocate(); diff --git a/src/CG-DNA/pair_oxdna_excv.cpp b/src/CG-DNA/pair_oxdna_excv.cpp index 2ae5995ebb..ae5752028a 100644 --- a/src/CG-DNA/pair_oxdna_excv.cpp +++ b/src/CG-DNA/pair_oxdna_excv.cpp @@ -39,7 +39,7 @@ PairOxdnaExcv::PairOxdnaExcv(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; writedata = 1; - + // set comm size needed by this Pair comm_forward = 9; } @@ -49,7 +49,7 @@ PairOxdnaExcv::PairOxdnaExcv(LAMMPS *lmp) : Pair(lmp) PairOxdnaExcv::~PairOxdnaExcv() { if (allocated) { - + memory->destroy(nx); memory->destroy(ny); memory->destroy(nz); @@ -127,7 +127,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag) // Cartesian unit vectors in lab frame double ax[3],ay[3],az[3]; double bx[3],by[3],bz[3]; - + double *special_lj = force->special_lj; double **x = atom->x; @@ -152,12 +152,12 @@ void PairOxdnaExcv::compute(int eflag, int vflag) alist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - - // loop over all local atoms, calculation of local reference frame + + // loop over all local atoms, calculation of local reference frame for (in = 0; in < atom->nlocal; in++) { - + int n = alist[in]; - double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame + double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame qn=bonus[ellipsoid[n]].quat; MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp); @@ -171,9 +171,9 @@ void PairOxdnaExcv::compute(int eflag, int vflag) nz[n][0] = nz_temp[0]; nz[n][1] = nz_temp[1]; nz[n][2] = nz_temp[2]; - + } - + comm->forward_comm_pair(this); // loop over pair interaction neighbors of my atoms @@ -181,7 +181,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag) for (ia = 0; ia < anum; ia++) { a = alist[ia]; - atype = type[a]; + atype = type[a]; ax[0] = nx[a][0]; ax[1] = nx[a][1]; @@ -191,7 +191,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag) ay[2] = ny[a][2]; az[0] = nz[a][0]; az[1] = nz[a][1]; - az[2] = nz[a][2]; + az[2] = nz[a][2]; // vector COM - backbone and base site a compute_interaction_sites(ax,ay,az,ra_cs,ra_cb); @@ -892,7 +892,7 @@ void PairOxdnaExcv::unpack_forward_comm(int n, int first, double *buf) nz[i][0] = buf[m++]; nz[i][1] = buf[m++]; nz[i][2] = buf[m++]; - } + } } /* ---------------------------------------------------------------------- */ @@ -900,7 +900,7 @@ void PairOxdnaExcv::unpack_forward_comm(int n, int first, double *buf) void *PairOxdnaExcv::extract(const char *str, int &dim) { dim = 2; - + if (strcmp(str,"nx") == 0) return (void *) nx; if (strcmp(str,"ny") == 0) return (void *) ny; if (strcmp(str,"nz") == 0) return (void *) nz; diff --git a/src/CG-DNA/pair_oxdna_hbond.cpp b/src/CG-DNA/pair_oxdna_hbond.cpp index eb04d74f2e..3d2e803aa8 100644 --- a/src/CG-DNA/pair_oxdna_hbond.cpp +++ b/src/CG-DNA/pair_oxdna_hbond.cpp @@ -178,8 +178,8 @@ void PairOxdnaHbond::compute(int eflag, int vflag) alist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - - // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv + + // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -195,7 +195,7 @@ void PairOxdnaHbond::compute(int eflag, int vflag) ax[0] = nx_xtrct[a][0]; ax[1] = nx_xtrct[a][1]; ax[2] = nx_xtrct[a][2]; - + ra_chb[0] = d_chb*ax[0]; ra_chb[1] = d_chb*ax[1]; ra_chb[2] = d_chb*ax[2]; @@ -214,7 +214,7 @@ void PairOxdnaHbond::compute(int eflag, int vflag) bx[0] = nx_xtrct[b][0]; bx[1] = nx_xtrct[b][1]; bx[2] = nx_xtrct[b][2]; - + rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; rb_chb[2] = d_chb*bx[2]; @@ -270,8 +270,8 @@ void PairOxdnaHbond::compute(int eflag, int vflag) b_hb3[atype][btype], dtheta_hb3_c[atype][btype]); // early rejection criterium - if (f4t3) { - + if (f4t3) { + az[0] = nz_xtrct[a][0]; az[1] = nz_xtrct[a][1]; az[2] = nz_xtrct[a][2]; diff --git a/src/CG-DNA/pair_oxdna_hbond.h b/src/CG-DNA/pair_oxdna_hbond.h index 1a138e82de..4afd9b8ad6 100644 --- a/src/CG-DNA/pair_oxdna_hbond.h +++ b/src/CG-DNA/pair_oxdna_hbond.h @@ -66,9 +66,9 @@ class PairOxdnaHbond : public Pair { double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast; double **b_hb8, **dtheta_hb8_c; - + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors - + int seqdepflag; virtual void allocate(); diff --git a/src/CG-DNA/pair_oxdna_stk.cpp b/src/CG-DNA/pair_oxdna_stk.cpp index 7665abbb79..feeb41d78b 100644 --- a/src/CG-DNA/pair_oxdna_stk.cpp +++ b/src/CG-DNA/pair_oxdna_stk.cpp @@ -254,8 +254,8 @@ void PairOxdnaStk::compute(int eflag, int vflag) evdwl = 0.0; ev_init(eflag,vflag); - - // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv + + // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -342,13 +342,13 @@ void PairOxdnaStk::compute(int eflag, int vflag) // early rejection criterium if (f1) { - + az[0] = nz_xtrct[a][0]; az[1] = nz_xtrct[a][1]; az[2] = nz_xtrct[a][2]; bz[0] = nz_xtrct[b][0]; bz[1] = nz_xtrct[b][1]; - bz[2] = nz_xtrct[b][2]; + bz[2] = nz_xtrct[b][2]; // theta4 angle and correction cost4 = MathExtra::dot3(bz,az); @@ -373,10 +373,10 @@ void PairOxdnaStk::compute(int eflag, int vflag) // early rejection criterium if (f4t5) { - + ay[0] = ny_xtrct[a][0]; ay[1] = ny_xtrct[a][1]; - ay[2] = ny_xtrct[a][2]; + ay[2] = ny_xtrct[a][2]; by[0] = ny_xtrct[b][0]; by[1] = ny_xtrct[b][1]; by[2] = ny_xtrct[b][2]; diff --git a/src/CG-DNA/pair_oxdna_stk.h b/src/CG-DNA/pair_oxdna_stk.h index 1f8a3a9249..5dbe9ea1e9 100644 --- a/src/CG-DNA/pair_oxdna_stk.h +++ b/src/CG-DNA/pair_oxdna_stk.h @@ -59,7 +59,7 @@ class PairOxdnaStk : public Pair { double **b_st6, **dtheta_st6_c; double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c; double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c; - + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors diff --git a/src/CG-DNA/pair_oxdna_xstk.cpp b/src/CG-DNA/pair_oxdna_xstk.cpp index bfb0d185a9..61dc1a8100 100644 --- a/src/CG-DNA/pair_oxdna_xstk.cpp +++ b/src/CG-DNA/pair_oxdna_xstk.cpp @@ -155,8 +155,8 @@ void PairOxdnaXstk::compute(int eflag, int vflag) alist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - - // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv + + // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); ny_xtrct = (double **) force->pair->extract("ny",dim); @@ -250,7 +250,7 @@ void PairOxdnaXstk::compute(int eflag, int vflag) // early rejection criterium if (f4t3) { - + az[0] = nz_xtrct[a][0]; az[1] = nz_xtrct[a][1]; az[2] = nz_xtrct[a][2]; diff --git a/src/CG-DNA/pair_oxdna_xstk.h b/src/CG-DNA/pair_oxdna_xstk.h index 2ad4c6027c..7695effc89 100644 --- a/src/CG-DNA/pair_oxdna_xstk.h +++ b/src/CG-DNA/pair_oxdna_xstk.h @@ -65,7 +65,7 @@ class PairOxdnaXstk : public Pair { double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast; double **b_xst8, **dtheta_xst8_c; - + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors diff --git a/src/CG-DNA/pair_oxrna2_stk.cpp b/src/CG-DNA/pair_oxrna2_stk.cpp index 18cd798fc7..c5fa47e25e 100644 --- a/src/CG-DNA/pair_oxrna2_stk.cpp +++ b/src/CG-DNA/pair_oxrna2_stk.cpp @@ -245,9 +245,9 @@ void PairOxrna2Stk::compute(int eflag, int vflag) double ra_cs[3],ra_cst[3]; double rb_cs[3],rb_cst[3]; - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; + // Cartesian unit vectors in lab frame + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -274,6 +274,12 @@ void PairOxrna2Stk::compute(int eflag, int vflag) evdwl = 0.0; ev_init(eflag,vflag); + // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + ny_xtrct = (double **) force->pair->extract("ny",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); + // loop over stacking interaction neighbors using bond topology for (in = 0; in < nbondlist; in++) { @@ -290,12 +296,26 @@ void PairOxrna2Stk::compute(int eflag, int vflag) } - // a now in 3' direction, b in 5' direction + // 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); + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; + ay[0] = ny_xtrct[a][0]; + ay[1] = ny_xtrct[a][1]; + ay[2] = ny_xtrct[a][2]; + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; + by[0] = ny_xtrct[b][0]; + by[1] = ny_xtrct[b][1]; + by[2] = ny_xtrct[b][2]; + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; // vector COM a - 5'-stacking site a ra_cst[0] = d_cst_x_5p*ax[0] + d_cst_y_5p*ay[0]; diff --git a/src/CG-DNA/pair_oxrna2_stk.h b/src/CG-DNA/pair_oxrna2_stk.h index 23650c98b5..e15614f126 100644 --- a/src/CG-DNA/pair_oxrna2_stk.h +++ b/src/CG-DNA/pair_oxrna2_stk.h @@ -61,6 +61,9 @@ class PairOxrna2Stk : public Pair { double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c; double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c; + // per-atom arrays for local unit vectors + double **nx_xtrct, **ny_xtrct, **nz_xtrct; + int seqdepflag; virtual void allocate(); diff --git a/src/CG-DNA/pair_oxrna2_xstk.cpp b/src/CG-DNA/pair_oxrna2_xstk.cpp index 8b44a9bde1..1f069e6b2f 100644 --- a/src/CG-DNA/pair_oxrna2_xstk.cpp +++ b/src/CG-DNA/pair_oxrna2_xstk.cpp @@ -120,9 +120,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag) // vectors COM-h-bonding site in lab frame double ra_chb[3],rb_chb[3]; - // quaternions and Cartesian unit vectors in lab frame - double *qa,ax[3],ay[3],az[3]; - double *qb,bx[3],by[3],bz[3]; + // Cartesian unit vectors in lab frame + double ax[3],ay[3],az[3]; + double bx[3],by[3],bz[3]; double **x = atom->x; double **f = atom->f; @@ -151,6 +151,11 @@ void PairOxrna2Xstk::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // n(x/z)_xtrct = extracted local unit vectors from oxdna_excv + int dim; + nx_xtrct = (double **) force->pair->extract("nx",dim); + nz_xtrct = (double **) force->pair->extract("nz",dim); + // loop over pair interaction neighbors of my atoms for (ia = 0; ia < anum; ia++) { @@ -158,8 +163,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag) a = alist[ia]; atype = type[a]; - qa=bonus[ellipsoid[a]].quat; - MathExtra::q_to_exyz(qa,ax,ay,az); + ax[0] = nx_xtrct[a][0]; + ax[1] = nx_xtrct[a][1]; + ax[2] = nx_xtrct[a][2]; ra_chb[0] = d_chb*ax[0]; ra_chb[1] = d_chb*ax[1]; @@ -176,8 +182,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag) btype = type[b]; - qb=bonus[ellipsoid[b]].quat; - MathExtra::q_to_exyz(qb,bx,by,bz); + bx[0] = nx_xtrct[b][0]; + bx[1] = nx_xtrct[b][1]; + bx[2] = nx_xtrct[b][2]; rb_chb[0] = d_chb*bx[0]; rb_chb[1] = d_chb*bx[1]; @@ -236,6 +243,10 @@ void PairOxrna2Xstk::compute(int eflag, int vflag) // early rejection criterium if (f4t3) { + az[0] = nz_xtrct[a][0]; + az[1] = nz_xtrct[a][1]; + az[2] = nz_xtrct[a][2]; + cost7 = -1.0*MathExtra::dot3(az,delr_hb_norm); if (cost7 > 1.0) cost7 = 1.0; if (cost7 < -1.0) cost7 = -1.0; @@ -250,6 +261,10 @@ void PairOxrna2Xstk::compute(int eflag, int vflag) // early rejection criterium if (f4t7) { + bz[0] = nz_xtrct[b][0]; + bz[1] = nz_xtrct[b][1]; + bz[2] = nz_xtrct[b][2]; + cost8 = MathExtra::dot3(bz,delr_hb_norm); if (cost8 > 1.0) cost8 = 1.0; if (cost8 < -1.0) cost8 = -1.0; diff --git a/src/CG-DNA/pair_oxrna2_xstk.h b/src/CG-DNA/pair_oxrna2_xstk.h index 187d0b67ae..6ece8e3f32 100644 --- a/src/CG-DNA/pair_oxrna2_xstk.h +++ b/src/CG-DNA/pair_oxrna2_xstk.h @@ -62,6 +62,9 @@ class PairOxrna2Xstk : public Pair { double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast; double **b_xst8, **dtheta_xst8_c; + // per-atom arrays for local unit vectors + double **nx_xtrct, **nz_xtrct; + virtual void allocate(); }; From d178a5ffa39a35ce67c2faa7c3e4b23826325e42 Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Wed, 17 Nov 2021 13:12:44 +0000 Subject: [PATCH 08/13] Extended test script --- examples/PACKAGES/cgdna/examples/test.sh | 224 +++++++++++++++++------ 1 file changed, 168 insertions(+), 56 deletions(-) diff --git a/examples/PACKAGES/cgdna/examples/test.sh b/examples/PACKAGES/cgdna/examples/test.sh index 795be9ef88..2e68ccab08 100755 --- a/examples/PACKAGES/cgdna/examples/test.sh +++ b/examples/PACKAGES/cgdna/examples/test.sh @@ -1,8 +1,9 @@ #! /bin/bash -DATE='2Jul21' -LMPDIR=/Users/ohenrich/Work/code/lammps +DATE='27Oct21' +TOL=1e-8 +LMPDIR=/Users/ohenrich/Work/code/lammps SRCDIR=$LMPDIR/src EXDIR=$LMPDIR/examples/PACKAGES/cgdna/examples @@ -11,6 +12,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then cd $SRCDIR make clean-all + make purge + make pu + make ps make -j8 mpi ###################################################### @@ -24,15 +28,27 @@ if [ $# -eq 1 ] && [ $1 = run ]; then mpirun -np 1 ./lmp_mpi < in.duplex1 > /dev/null mv log.lammps log.$DATE.duplex1.g++.1 - grep etot log.$DATE.duplex1.g++.1 > e_test.dat - grep etot ../log*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex1.g++.1 > e_test.1.dat + grep etot ../log*1 > e_old.1.dat + ndiff -relerr $TOL e_test.1.dat e_old.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task passed" + else + echo "# 1 MPI-task unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.duplex1 > /dev/null mv log.lammps log.$DATE.duplex1.g++.4 - grep etot log.$DATE.duplex1.g++.4 > e_test.dat - grep etot ../log*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex1.g++.4 > e_test.4.dat + grep etot ../log*4 > e_old.4.dat + ndiff -relerr $TOL e_test.4.dat e_old.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks passed" + else + echo "# 4 MPI-tasks unsuccessful" + fi ###################################################### @@ -47,15 +63,27 @@ if [ $# -eq 1 ] && [ $1 = run ]; then mpirun -np 1 ./lmp_mpi < in.duplex2 > /dev/null mv log.lammps log.$DATE.duplex2.g++.1 - grep etot log.$DATE.duplex2.g++.1 > e_test.dat - grep etot ../log*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex2.g++.1 > e_test.1.dat + grep etot ../log*1 > e_old.1.dat + ndiff -relerr $TOL e_test.1.dat e_old.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task passed" + else + echo "# 1 MPI-task unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null mv log.lammps log.$DATE.duplex2.g++.4 - grep etot log.$DATE.duplex2.g++.4 > e_test.dat - grep etot ../log*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex2.g++.4 > e_test.4.dat + grep etot ../log*4 > e_old.4.dat + ndiff -relerr $TOL e_test.4.dat e_old.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks passed" + else + echo "# 4 MPI-tasks unsuccessful" + fi ###################################################### @@ -70,15 +98,27 @@ if [ $# -eq 1 ] && [ $1 = run ]; then mpirun -np 1 ./lmp_mpi < in.duplex1 > /dev/null mv log.lammps log.$DATE.duplex1.g++.1 - grep etot log.$DATE.duplex1.g++.1 > e_test.dat - grep etot ../log*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex1.g++.1 > e_test.1.dat + grep etot ../log*1 > e_old.1.dat + ndiff -relerr $TOL e_test.1.dat e_old.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task passed" + else + echo "# 1 MPI-task unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.duplex1 > /dev/null mv log.lammps log.$DATE.duplex1.g++.4 - grep etot log.$DATE.duplex1.g++.4 > e_test.dat - grep etot ../log*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex1.g++.4 > e_test.4.dat + grep etot ../log*4 > e_old.4.dat + ndiff -relerr $TOL e_test.4.dat e_old.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks passed" + else + echo "# 4 MPI-tasks unsuccessful" + fi ###################################################### @@ -93,15 +133,27 @@ if [ $# -eq 1 ] && [ $1 = run ]; then mpirun -np 1 ./lmp_mpi < in.duplex2 > /dev/null mv log.lammps log.$DATE.duplex2.g++.1 - grep etot log.$DATE.duplex2.g++.1 > e_test.dat - grep etot ../log*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex2.g++.1 > e_test.1.dat + grep etot ../log*1 > e_old.1.dat + ndiff -relerr $TOL e_test.1.dat e_old.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task passed" + else + echo "# 1 MPI-task unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null mv log.lammps log.$DATE.duplex2.g++.4 - grep etot log.$DATE.duplex2.g++.4 > e_test.dat - grep etot ../log*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex2.g++.4 > e_test.4.dat + grep etot ../log*4 > e_old.4.dat + ndiff -relerr $TOL e_test.4.dat e_old.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks passed" + else + echo "# 4 MPI-tasks unsuccessful" + fi ###################################################### @@ -116,15 +168,27 @@ if [ $# -eq 1 ] && [ $1 = run ]; then mpirun -np 1 ./lmp_mpi < in.duplex3 > /dev/null mv log.lammps log.$DATE.duplex3.g++.1 - grep etot log.$DATE.duplex3.g++.1 > e_test.dat - grep etot ../log*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex3.g++.1 > e_test.1.dat + grep etot ../log*1 > e_old.1.dat + ndiff -relerr $TOL e_test.1.dat e_old.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task passed" + else + echo "# 1 MPI-task unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.duplex3 > /dev/null mv log.lammps log.$DATE.duplex3.g++.4 - grep etot log.$DATE.duplex3.g++.4 > e_test.dat - grep etot ../log*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex3.g++.4 > e_test.4.dat + grep etot ../log*4 > e_old.4.dat + ndiff -relerr $TOL e_test.4.dat e_old.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks passed" + else + echo "# 4 MPI-tasks unsuccessful" + fi ###################################################### @@ -141,27 +205,51 @@ if [ $# -eq 1 ] && [ $1 = run ]; then mpirun -np 1 ./lmp_mpi < in.duplex4.4type > /dev/null mv log.lammps log.$DATE.duplex4.4type.g++.1 - grep etot log.$DATE.duplex4.4type.g++.1 > e_test.dat - grep etot ../log*4type*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex4.4type.g++.1 > e_test.4type.1.dat + grep etot ../log*4type*1 > e_old.4type.1.dat + ndiff -relerr $TOL e_test.4type.1.dat e_old.4type.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task 4 types passed" + else + echo "# 1 MPI-task 4 types unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.duplex4.4type > /dev/null mv log.lammps log.$DATE.duplex4.4type.g++.4 - grep etot log.$DATE.duplex4.4type.g++.4 > e_test.dat - grep etot ../log*4type*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex4.4type.g++.4 > e_test.4type.4.dat + grep etot ../log*4type*4 > e_old.4type.4.dat + ndiff -relerr $TOL e_test.4type.4.dat e_old.4type.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks 4 types passed" + else + echo "# 4 MPI-tasks 4 types unsuccessful" + fi mpirun -np 1 ./lmp_mpi < in.duplex4.8type > /dev/null mv log.lammps log.$DATE.duplex4.8type.g++.1 - grep etot log.$DATE.duplex4.8type.g++.1 > e_test.dat - grep etot ../log*8type*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex4.8type.g++.1 > e_test.8type.1.dat + grep etot ../log*8type*1 > e_old.8type.1.dat + ndiff -relerr $TOL e_test.8type.1.dat e_old.8type.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task 8 types passed" + else + echo "# 1 MPI-task 8 types unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.duplex4.8type > /dev/null mv log.lammps log.$DATE.duplex4.8type.g++.4 - grep etot log.$DATE.duplex4.8type.g++.4 > e_test.dat - grep etot ../log*8type*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex4.8type.g++.4 > e_test.8type.4.dat + grep etot ../log*8type*4 > e_old.8type.4.dat + ndiff -relerr $TOL e_test.8type.4.dat e_old.8type.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks 8 types passed" + else + echo "# 4 MPI-tasks 8 types unsuccessful" + fi ###################################################### @@ -176,15 +264,27 @@ if [ $# -eq 1 ] && [ $1 = run ]; then mpirun -np 1 ./lmp_mpi < in.dsring > /dev/null mv log.lammps log.$DATE.dsring.g++.1 - grep etot log.$DATE.dsring.g++.1 > e_test.dat - grep etot ../log*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.dsring.g++.1 > e_test.1.dat + grep etot ../log*1 > e_old.1.dat + ndiff -relerr $TOL e_test.1.dat e_old.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task passed" + else + echo "# 1 MPI-task unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.dsring > /dev/null mv log.lammps log.$DATE.dsring.g++.4 - grep etot log.$DATE.dsring.g++.4 > e_test.dat - grep etot ../log*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.dsring.g++.4 > e_test.4.dat + grep etot ../log*4 > e_old.4.dat + ndiff -relerr $TOL e_test.4.dat e_old.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks passed" + else + echo "# 4 MPI-tasks unsuccessful" + fi ###################################################### @@ -199,15 +299,27 @@ if [ $# -eq 1 ] && [ $1 = run ]; then mpirun -np 1 ./lmp_mpi < in.duplex2 > /dev/null mv log.lammps log.$DATE.duplex2.g++.1 - grep etot log.$DATE.duplex2.g++.1 > e_test.dat - grep etot ../log*1 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex2.g++.1 > e_test.1.dat + grep etot ../log*1 > e_old.1.dat + ndiff -relerr $TOL e_test.1.dat e_old.1.dat + if [ $? -eq 0 ]; + then + echo "# 1 MPI-task passed" + else + echo "# 1 MPI-task unsuccessful" + fi mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null mv log.lammps log.$DATE.duplex2.g++.4 - grep etot log.$DATE.duplex2.g++.4 > e_test.dat - grep etot ../log*4 > e_old.dat - ndiff -relerr 1e-8 e_test.dat e_old.dat + grep etot log.$DATE.duplex2.g++.4 > e_test.4.dat + grep etot ../log*4 > e_old.4.dat + ndiff -relerr $TOL e_test.4.dat e_old.4.dat + if [ $? -eq 0 ]; + then + echo "# 4 MPI-tasks passed" + else + echo "# 4 MPI-tasks unsuccessful" + fi ###################################################### echo '# Done' From 1f81e2afad621c426472bd5753d4ead30121d70d Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Fri, 7 Jan 2022 13:45:54 +0000 Subject: [PATCH 09/13] Added duplication of stdout into logfile --- examples/PACKAGES/cgdna/examples/test.sh | 105 ++++++++++++----------- 1 file changed, 53 insertions(+), 52 deletions(-) diff --git a/examples/PACKAGES/cgdna/examples/test.sh b/examples/PACKAGES/cgdna/examples/test.sh index 2e68ccab08..152047b94b 100755 --- a/examples/PACKAGES/cgdna/examples/test.sh +++ b/examples/PACKAGES/cgdna/examples/test.sh @@ -1,6 +1,6 @@ #! /bin/bash -DATE='27Oct21' +DATE='14Dec21' TOL=1e-8 LMPDIR=/Users/ohenrich/Work/code/lammps @@ -8,17 +8,17 @@ SRCDIR=$LMPDIR/src EXDIR=$LMPDIR/examples/PACKAGES/cgdna/examples if [ $# -eq 1 ] && [ $1 = run ]; then - echo '# Compiling executable in' $SRCDIR + echo '# Compiling executable in' $SRCDIR | tee -a $EXDIR/test.log cd $SRCDIR - make clean-all - make purge - make pu - make ps - make -j8 mpi + make clean-all | tee -a $EXDIR/test.log + make purge | tee -a $EXDIR/test.log + make pu | tee -a $EXDIR/test.log + make ps | tee -a $EXDIR/test.log + make -j8 mpi | tee -a $EXDIR/test.log ###################################################### - echo '# Running oxDNA duplex1 test' + echo '# Running oxDNA duplex1 test' | tee -a $EXDIR/test.log cd $EXDIR/oxDNA/duplex1 mkdir test cd test @@ -33,9 +33,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.1.dat e_old.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task passed" + echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task unsuccessful" + echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.duplex1 > /dev/null @@ -45,15 +45,15 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4.dat e_old.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks passed" + echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks unsuccessful" + echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log fi ###################################################### ###################################################### - echo '# Running oxDNA duplex2 test' + echo '# Running oxDNA duplex2 test' | tee -a $EXDIR/test.log cd $EXDIR/oxDNA/duplex2 mkdir test cd test @@ -68,9 +68,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.1.dat e_old.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task passed" + echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task unsuccessful" + echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null @@ -80,15 +80,15 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4.dat e_old.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks passed" + echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks unsuccessful" + echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log fi ###################################################### ###################################################### - echo '# Running oxDNA2 duplex1 test' + echo '# Running oxDNA2 duplex1 test' | tee -a $EXDIR/test.log cd $EXDIR/oxDNA2/duplex1 mkdir test cd test @@ -103,9 +103,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.1.dat e_old.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task passed" + echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task unsuccessful" + echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.duplex1 > /dev/null @@ -115,15 +115,15 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4.dat e_old.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks passed" + echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks unsuccessful" + echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log fi ###################################################### ###################################################### - echo '# Running oxDNA2 duplex2 test' + echo '# Running oxDNA2 duplex2 test' | tee -a $EXDIR/test.log cd $EXDIR/oxDNA2/duplex2 mkdir test cd test @@ -138,9 +138,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.1.dat e_old.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task passed" + echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task unsuccessful" + echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null @@ -150,15 +150,15 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4.dat e_old.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks passed" + echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks unsuccessful" + echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log fi ###################################################### ###################################################### - echo '# Running oxDNA2 duplex3 test' + echo '# Running oxDNA2 duplex3 test' | tee -a $EXDIR/test.log cd $EXDIR/oxDNA2/duplex3 mkdir test cd test @@ -173,9 +173,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.1.dat e_old.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task passed" + echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task unsuccessful" + echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.duplex3 > /dev/null @@ -185,15 +185,15 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4.dat e_old.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks passed" + echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks unsuccessful" + echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log fi ###################################################### ###################################################### - echo '# Running oxDNA2 unique_bp test' + echo '# Running oxDNA2 unique_bp test' | tee -a $EXDIR/test.log cd $EXDIR/oxDNA2/unique_bp mkdir test cd test @@ -210,9 +210,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4type.1.dat e_old.4type.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task 4 types passed" + echo "# 1 MPI-task 4 types passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task 4 types unsuccessful" + echo "# 1 MPI-task 4 types unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.duplex4.4type > /dev/null @@ -222,9 +222,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4type.4.dat e_old.4type.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks 4 types passed" + echo "# 4 MPI-tasks 4 types passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks 4 types unsuccessful" + echo "# 4 MPI-tasks 4 types unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 1 ./lmp_mpi < in.duplex4.8type > /dev/null @@ -234,9 +234,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.8type.1.dat e_old.8type.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task 8 types passed" + echo "# 1 MPI-task 8 types passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task 8 types unsuccessful" + echo "# 1 MPI-task 8 types unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.duplex4.8type > /dev/null @@ -246,15 +246,15 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.8type.4.dat e_old.8type.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks 8 types passed" + echo "# 4 MPI-tasks 8 types passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks 8 types unsuccessful" + echo "# 4 MPI-tasks 8 types unsuccessful" | tee -a $EXDIR/test.log fi ###################################################### ###################################################### - echo '# Running oxDNA2 dsring test' + echo '# Running oxDNA2 dsring test' | tee -a $EXDIR/test.log cd $EXDIR/oxDNA2/dsring mkdir test cd test @@ -269,9 +269,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.1.dat e_old.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task passed" + echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task unsuccessful" + echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.dsring > /dev/null @@ -281,15 +281,15 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4.dat e_old.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks passed" + echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks unsuccessful" + echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log fi ###################################################### ###################################################### - echo '# Running oxRNA2 duplex2 test' + echo '# Running oxRNA2 duplex2 test' | tee -a $EXDIR/test.log cd $EXDIR/oxRNA2/duplex2 mkdir test cd test @@ -304,9 +304,9 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.1.dat e_old.1.dat if [ $? -eq 0 ]; then - echo "# 1 MPI-task passed" + echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log else - echo "# 1 MPI-task unsuccessful" + echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log fi mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null @@ -316,13 +316,13 @@ if [ $# -eq 1 ] && [ $1 = run ]; then ndiff -relerr $TOL e_test.4.dat e_old.4.dat if [ $? -eq 0 ]; then - echo "# 4 MPI-tasks passed" + echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log else - echo "# 4 MPI-tasks unsuccessful" + echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log fi ###################################################### - echo '# Done' + echo '# Done' | tee -a $EXDIR/test.log elif [ $# -eq 1 ] && [ $1 = clean ]; then echo '# Deleting test directories' @@ -334,6 +334,7 @@ elif [ $# -eq 1 ] && [ $1 = clean ]; then rm -rf $EXDIR/oxDNA2/unique_bp/test rm -rf $EXDIR/oxDNA2/dsring/test rm -rf $EXDIR/oxRNA2/duplex2/test + rm -rf $EXDIR/test.log echo '# Done' else From 9b6bb00d07e938ce69e683734a74e52dfe620a01 Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Wed, 26 Jan 2022 19:20:34 +0000 Subject: [PATCH 10/13] Added unit test for atom_style oxdna --- unittest/formats/test_atom_styles.cpp | 444 ++++++++++++++++++++++++++ 1 file changed, 444 insertions(+) diff --git a/unittest/formats/test_atom_styles.cpp b/unittest/formats/test_atom_styles.cpp index 59c0e1350c..5076e92ca3 100644 --- a/unittest/formats/test_atom_styles.cpp +++ b/unittest/formats/test_atom_styles.cpp @@ -4809,6 +4809,450 @@ TEST_F(AtomStyleTest, property_atom) EXPECT_NEAR(three[GETIDX(2)], 0.5, EPSILON); } +TEST_F(AtomStyleTest, oxdna) +{ + if (!LAMMPS::is_installed_pkg("MOLECULE")) GTEST_SKIP(); + if (!LAMMPS::is_installed_pkg("ASPHERE")) GTEST_SKIP(); + if (!LAMMPS::is_installed_pkg("CG-DNA")) GTEST_SKIP(); + + BEGIN_HIDE_OUTPUT(); + command("atom_style hybrid bond ellipsoid oxdna"); + END_HIDE_OUTPUT(); + + AtomState expected; + expected.atom_style = "hybrid"; + expected.molecular = Atom::MOLECULAR; + expected.tag_enable = 1; + expected.molecule_flag = 1; + expected.ellipsoid_flag = 1; + expected.rmass_flag = 1; + expected.torque_flag = 1; + expected.angmom_flag = 1; + expected.has_type = true; + expected.has_mask = true; + expected.has_image = true; + expected.has_x = true; + expected.has_v = true; + expected.has_f = true; + expected.has_bonds = true; + expected.has_nspecial = true; + expected.has_special = true; + expected.map_style = 3; + + ASSERT_ATOM_STATE_EQ(lmp->atom, expected); + + auto hybrid = (AtomVecHybrid *)lmp->atom->avec; + ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid")); + ASSERT_EQ(hybrid->nstyles, 3); + ASSERT_THAT(std::string(hybrid->keywords[0]), Eq("bond")); + ASSERT_THAT(std::string(hybrid->keywords[1]), Eq("ellipsoid")); + ASSERT_THAT(std::string(hybrid->keywords[2]), Eq("oxdna")); + ASSERT_NE(hybrid->styles[0], nullptr); + ASSERT_NE(hybrid->styles[1], nullptr); + ASSERT_NE(hybrid->styles[2], nullptr); + + BEGIN_HIDE_OUTPUT(); + command("units lj"); + command("dimension 3"); + command("newton on"); + command("boundary p p p"); + command("atom_modify sort 0 1.0"); + command("neighbor 2.0 bin"); + command("neigh_modify every 1 delay 0 check yes"); + command("region mybox block -20 20 -20 20 -20 20"); + command("create_box 4 mybox bond/types 1 extra/bond/per/atom 2 extra/special/per/atom 4"); + command("create_atoms 1 single -0.33741452300167507 -0.43708835412476305 0.6450685042019271"); + command("create_atoms 2 single -0.32142606102826937 -0.7137743037592722 1.1817366147004618"); + command("create_atoms 3 single -0.130363628207774 -0.9147144801536078 1.62581312195109"); + command("create_atoms 4 single 0.16795127962282844 -0.9808507459807022 2.0894908590909003"); + command("create_atoms 1 single 0.46370423490634166 -0.7803347954883079 2.4251986815515827"); + command("create_atoms 4 single -0.4462950185476711 0.09062163051035639 2.4668941268777607"); + command("create_atoms 1 single -0.03377054097560965 0.20979847489755046 2.078208732038921"); + command("create_atoms 2 single 0.3297325391466579 0.17657587120899895 1.7206328374934152"); + command("create_atoms 3 single 0.6063699309305985 0.04682595158675571 1.2335049647817748"); + command("create_atoms 4 single 0.8003979559814726 -0.364393011459011 0.9884025318908612"); + command("set type 1 mass 3.1575"); + command("set type 2 mass 3.1575"); + command("set type 3 mass 3.1575"); + command("set type 4 mass 3.1575"); + command("mass 1 3.1575"); + command("mass 2 3.1575"); + command("mass 3 3.1575"); + command("mass 4 3.1575"); + command("set atom 1 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 2 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 3 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 4 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 5 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 6 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 7 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 8 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 9 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 10 shape 1.173984503142341 1.173984503142341 1.173984503142341"); + command("set atom 1 quat 0.120438343611269 -0.970540441176996 0.208676441957758 16.990727782866998"); + command("set atom 2 quat 0.122039415796829 -0.068232256412985 0.990177125658213 40.001729435287870"); + command("set atom 3 quat 0.052760168698289 0.030943512185297 0.998127679033382 69.627682451632380"); + command("set atom 4 quat -0.037622918613871 0.030623545471522 0.998822664169035 97.038820280300570"); + command("set atom 5 quat 0.055056042946138 0.077631917807377 0.995460756369964 137.7813218321917"); + command("set atom 6 quat 0.931128471673637 -0.355724722922553 -0.080372201291206 166.2836226291888"); + command("set atom 7 quat 0.753526078198930 -0.648440397941919 0.108275111595674 200.6802564250672"); + command("set atom 8 quat 0.553942138074214 -0.829580511279186 0.070315595507185 192.0355407659524"); + command("set atom 9 quat -0.373540155765431 0.913070802138105 -0.163613759548524 171.0789308260751"); + command("set atom 10 quat 0.027515673832457 0.998248649922676 -0.052369080773879 161.2621224558284"); + command("bond_style oxdna2/fene"); + command("bond_coeff * 2.0 0.25 0.7564"); + command("special_bonds lj 0 1 1"); + command("create_bonds single/bond 1 1 2"); + command("create_bonds single/bond 1 2 3"); + command("create_bonds single/bond 1 3 4"); + command("create_bonds single/bond 1 4 5"); + command("create_bonds single/bond 1 6 7"); + command("create_bonds single/bond 1 7 8"); + command("create_bonds single/bond 1 8 9"); + command("create_bonds single/bond 1 9 10"); + command("pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh"); + command("pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32"); + command("pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65"); + command("pair_coeff * * oxdna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"); + command("pair_coeff 1 4 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"); + command("pair_coeff 2 3 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"); + command("pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68"); + command("pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793"); + command("pair_coeff * * oxdna2/dh 0.1 0.2 0.815"); + END_HIDE_OUTPUT(); + + ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid")); + ASSERT_NE(lmp->atom->avec, nullptr); + ASSERT_EQ(lmp->atom->natoms, 10); + ASSERT_EQ(lmp->atom->nbonds, 8); + ASSERT_EQ(lmp->atom->nbondtypes, 1); + ASSERT_EQ(lmp->atom->nellipsoids, 10); + ASSERT_EQ(lmp->atom->nlocal, 10); + ASSERT_EQ(lmp->atom->nghost, 0); + ASSERT_NE(lmp->atom->nmax, -1); + ASSERT_EQ(lmp->atom->tag_enable, 1); + ASSERT_EQ(lmp->atom->molecular, Atom::MOLECULAR); + ASSERT_EQ(lmp->atom->ntypes, 4); + ASSERT_EQ(lmp->atom->nextra_grow, 0); + ASSERT_EQ(lmp->atom->nextra_restart, 0); + ASSERT_EQ(lmp->atom->nextra_border, 0); + ASSERT_EQ(lmp->atom->nextra_grow_max, 0); + ASSERT_EQ(lmp->atom->nextra_restart_max, 0); + ASSERT_EQ(lmp->atom->nextra_border_max, 0); + ASSERT_EQ(lmp->atom->nextra_store, 0); + ASSERT_EQ(lmp->atom->extra_grow, nullptr); + ASSERT_EQ(lmp->atom->extra_restart, nullptr); + ASSERT_EQ(lmp->atom->extra_border, nullptr); + ASSERT_EQ(lmp->atom->extra, nullptr); + ASSERT_NE(lmp->atom->mass, nullptr); + ASSERT_NE(lmp->atom->rmass, nullptr); + ASSERT_EQ(lmp->atom->ellipsoid_flag, 1); + ASSERT_NE(lmp->atom->ellipsoid, nullptr); + ASSERT_NE(lmp->atom->mass_setflag, nullptr); + ASSERT_NE(lmp->atom->id5p, nullptr); + + BEGIN_HIDE_OUTPUT(); + command("write_data test_atom_styles.data nocoeff"); + command("clear"); + command("units lj"); + command("dimension 3"); + command("newton on"); + command("boundary p p p"); + command("atom_style hybrid bond ellipsoid oxdna"); + command("atom_modify sort 0 1.0"); + command("neighbor 2.0 bin"); + command("neigh_modify every 1 delay 0 check yes"); + command("read_data test_atom_styles.data"); + command("set type 1 mass 3.1575"); + command("set type 2 mass 3.1575"); + command("set type 3 mass 3.1575"); + command("set type 4 mass 3.1575"); + command("mass 1 3.1575"); + command("mass 2 3.1575"); + command("mass 3 3.1575"); + command("mass 4 3.1575"); + command("bond_style oxdna2/fene"); + command("bond_coeff * 2.0 0.25 0.7564"); + command("special_bonds lj 0 1 1"); + command("pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh"); + command("pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32"); + command("pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65"); + command("pair_coeff * * oxdna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"); + command("pair_coeff 1 4 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"); + command("pair_coeff 2 3 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"); + command("pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68"); + command("pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793"); + command("pair_coeff * * oxdna2/dh 0.1 0.2 0.815"); + + ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid")); + ASSERT_NE(lmp->atom->avec, nullptr); + hybrid = (AtomVecHybrid *)lmp->atom->avec; + + ASSERT_EQ(hybrid->nstyles, 3); + ASSERT_THAT(std::string(hybrid->keywords[0]), Eq("bond")); + ASSERT_THAT(std::string(hybrid->keywords[1]), Eq("ellipsoid")); + ASSERT_THAT(std::string(hybrid->keywords[2]), Eq("oxdna")); + ASSERT_NE(hybrid->styles[0], nullptr); + ASSERT_NE(hybrid->styles[1], nullptr); + ASSERT_NE(hybrid->styles[2], nullptr); + + ASSERT_EQ(lmp->atom->natoms, 10); + ASSERT_EQ(lmp->atom->nbonds, 8); + ASSERT_EQ(lmp->atom->nbondtypes, 1); + ASSERT_EQ(lmp->atom->nellipsoids, 10); + ASSERT_EQ(lmp->atom->nlocal, 10); + ASSERT_EQ(lmp->atom->nghost, 0); + ASSERT_NE(lmp->atom->nmax, -1); + ASSERT_EQ(lmp->atom->tag_enable, 1); + ASSERT_EQ(lmp->atom->molecular, Atom::MOLECULAR); + ASSERT_EQ(lmp->atom->ntypes, 4); + ASSERT_EQ(lmp->atom->nextra_grow, 0); + ASSERT_EQ(lmp->atom->nextra_restart, 0); + ASSERT_EQ(lmp->atom->nextra_border, 0); + ASSERT_EQ(lmp->atom->nextra_grow_max, 0); + ASSERT_EQ(lmp->atom->nextra_restart_max, 0); + ASSERT_EQ(lmp->atom->nextra_border_max, 0); + ASSERT_EQ(lmp->atom->nextra_store, 0); + ASSERT_EQ(lmp->atom->extra_grow, nullptr); + ASSERT_EQ(lmp->atom->extra_restart, nullptr); + ASSERT_EQ(lmp->atom->extra_border, nullptr); + ASSERT_EQ(lmp->atom->extra, nullptr); + ASSERT_NE(lmp->atom->mass, nullptr); + ASSERT_NE(lmp->atom->rmass, nullptr); + ASSERT_EQ(lmp->atom->ellipsoid_flag, 1); + ASSERT_NE(lmp->atom->ellipsoid, nullptr); + ASSERT_NE(lmp->atom->mass_setflag, nullptr); + ASSERT_NE(lmp->atom->id5p, nullptr); + + auto x = lmp->atom->x; + auto v = lmp->atom->v; + auto type = lmp->atom->type; + auto ellipsoid = lmp->atom->ellipsoid; + auto rmass = lmp->atom->rmass; + + auto avec = (AtomVecEllipsoid *)hybrid->styles[1]; + auto bonus = avec->bonus; + + EXPECT_NEAR(x[GETIDX(1)][0], -0.33741452300167507, EPSILON); + EXPECT_NEAR(x[GETIDX(1)][1], -0.43708835412476305, EPSILON); + EXPECT_NEAR(x[GETIDX(1)][2], 0.6450685042019271, EPSILON); + EXPECT_NEAR(x[GETIDX(2)][0], -0.32142606102826937, EPSILON); + EXPECT_NEAR(x[GETIDX(2)][1], -0.7137743037592722, EPSILON); + EXPECT_NEAR(x[GETIDX(2)][2], 1.1817366147004618, EPSILON); + EXPECT_NEAR(x[GETIDX(3)][0], -0.130363628207774, EPSILON); + EXPECT_NEAR(x[GETIDX(3)][1], -0.9147144801536078, EPSILON); + EXPECT_NEAR(x[GETIDX(3)][2], 1.62581312195109, EPSILON); + EXPECT_NEAR(x[GETIDX(4)][0], 0.16795127962282844, EPSILON); + EXPECT_NEAR(x[GETIDX(4)][1], -0.9808507459807022, EPSILON); + EXPECT_NEAR(x[GETIDX(4)][2], 2.0894908590909003, EPSILON); + EXPECT_NEAR(x[GETIDX(5)][0], 0.46370423490634166, EPSILON); + EXPECT_NEAR(x[GETIDX(5)][1], -0.7803347954883079, EPSILON); + EXPECT_NEAR(x[GETIDX(5)][2], 2.4251986815515827, EPSILON); + EXPECT_NEAR(x[GETIDX(6)][0], -0.4462950185476711, EPSILON); + EXPECT_NEAR(x[GETIDX(6)][1], 0.09062163051035639, EPSILON); + EXPECT_NEAR(x[GETIDX(6)][2], 2.4668941268777607, EPSILON); + EXPECT_NEAR(x[GETIDX(7)][0], -0.03377054097560965, EPSILON); + EXPECT_NEAR(x[GETIDX(7)][1], 0.20979847489755046, EPSILON); + EXPECT_NEAR(x[GETIDX(7)][2], 2.078208732038921, EPSILON); + EXPECT_NEAR(x[GETIDX(8)][0], 0.3297325391466579, EPSILON); + EXPECT_NEAR(x[GETIDX(8)][1], 0.17657587120899895, EPSILON); + EXPECT_NEAR(x[GETIDX(8)][2], 1.7206328374934152, EPSILON); + EXPECT_NEAR(x[GETIDX(9)][0], 0.6063699309305985, EPSILON); + EXPECT_NEAR(x[GETIDX(9)][1], 0.04682595158675571, EPSILON); + EXPECT_NEAR(x[GETIDX(9)][2], 1.2335049647817748, EPSILON); + EXPECT_NEAR(x[GETIDX(10)][0], 0.8003979559814726, EPSILON); + EXPECT_NEAR(x[GETIDX(10)][1], -0.364393011459011, EPSILON); + EXPECT_NEAR(x[GETIDX(10)][2], 0.9884025318908612, EPSILON); + + EXPECT_NEAR(v[GETIDX(1)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(1)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(1)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(2)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(2)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(2)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(3)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(3)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(3)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(4)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(4)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(4)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(5)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(5)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(5)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(6)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(6)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(6)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(7)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(7)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(7)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(8)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(8)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(8)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(9)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(9)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(9)][2], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(10)][0], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(10)][1], 0.0, EPSILON); + EXPECT_NEAR(v[GETIDX(10)][2], 0.0, EPSILON); + + ASSERT_EQ(type[GETIDX(1)], 1); + ASSERT_EQ(type[GETIDX(2)], 2); + ASSERT_EQ(type[GETIDX(3)], 3); + ASSERT_EQ(type[GETIDX(4)], 4); + ASSERT_EQ(type[GETIDX(5)], 1); + ASSERT_EQ(type[GETIDX(6)], 4); + ASSERT_EQ(type[GETIDX(7)], 1); + ASSERT_EQ(type[GETIDX(8)], 2); + ASSERT_EQ(type[GETIDX(9)], 3); + ASSERT_EQ(type[GETIDX(10)], 4); + + ASSERT_EQ(ellipsoid[GETIDX(1)], 0); + ASSERT_EQ(ellipsoid[GETIDX(2)], 1); + ASSERT_EQ(ellipsoid[GETIDX(3)], 2); + ASSERT_EQ(ellipsoid[GETIDX(4)], 3); + ASSERT_EQ(ellipsoid[GETIDX(5)], 4); + ASSERT_EQ(ellipsoid[GETIDX(6)], 5); + ASSERT_EQ(ellipsoid[GETIDX(7)], 6); + ASSERT_EQ(ellipsoid[GETIDX(8)], 7); + ASSERT_EQ(ellipsoid[GETIDX(9)], 8); + ASSERT_EQ(ellipsoid[GETIDX(10)], 9); + + EXPECT_NEAR(rmass[GETIDX(1)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(2)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(3)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(4)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(5)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(6)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(7)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(8)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(9)], 3.1575, EPSILON); + EXPECT_NEAR(rmass[GETIDX(10)], 3.1575, EPSILON); + + EXPECT_NEAR(bonus[0].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[0].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[0].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[1].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[1].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[1].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[2].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[2].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[2].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[3].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[3].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[3].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[4].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[4].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[4].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[5].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[5].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[5].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[6].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[6].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[6].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[7].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[7].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[7].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[8].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[8].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[8].shape[2], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[9].shape[0], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[9].shape[1], 0.5869922515711705, EPSILON); + EXPECT_NEAR(bonus[9].shape[2], 0.5869922515711705, EPSILON); + + EXPECT_NEAR(bonus[0].quat[0], 0.9890278201757743, EPSILON); + EXPECT_NEAR(bonus[0].quat[1], 0.01779228232037064, EPSILON); + EXPECT_NEAR(bonus[0].quat[2], -0.14337734159225404, EPSILON); + EXPECT_NEAR(bonus[0].quat[3], 0.030827642240801516, EPSILON); + EXPECT_NEAR(bonus[1].quat[0], 0.939687458852748, EPSILON); + EXPECT_NEAR(bonus[1].quat[1], 0.04174166924055095, EPSILON); + EXPECT_NEAR(bonus[1].quat[2], -0.023337773785056866, EPSILON); + EXPECT_NEAR(bonus[1].quat[3], 0.338674565089608, EPSILON); + EXPECT_NEAR(bonus[2].quat[0], 0.8210113150655425, EPSILON); + EXPECT_NEAR(bonus[2].quat[1], 0.03012140921736572, EPSILON); + EXPECT_NEAR(bonus[2].quat[2], 0.017666019956944813, EPSILON); + EXPECT_NEAR(bonus[2].quat[3], 0.5698429897612057, EPSILON); + EXPECT_NEAR(bonus[3].quat[0], 0.6623662858285051, EPSILON); + EXPECT_NEAR(bonus[3].quat[1], -0.028186343967346823, EPSILON); + EXPECT_NEAR(bonus[3].quat[2], 0.022942552517501488, EPSILON); + EXPECT_NEAR(bonus[3].quat[3], 0.7482981175276918, EPSILON); + EXPECT_NEAR(bonus[4].quat[0], 0.3601488726765216, EPSILON); + EXPECT_NEAR(bonus[4].quat[1], 0.0513614985821682, EPSILON); + EXPECT_NEAR(bonus[4].quat[2], 0.0724224158335286, EPSILON); + EXPECT_NEAR(bonus[4].quat[3], 0.9286602067807472, EPSILON); + EXPECT_NEAR(bonus[5].quat[0], 0.11941234710084649, EPSILON); + EXPECT_NEAR(bonus[5].quat[1], 0.9244660117493703, EPSILON); + EXPECT_NEAR(bonus[5].quat[2], -0.35317942248051865, EPSILON); + EXPECT_NEAR(bonus[5].quat[3], -0.07979711784524246, EPSILON); + EXPECT_NEAR(bonus[6].quat[0], -0.17949125421205164, EPSILON); + EXPECT_NEAR(bonus[6].quat[1], 0.7412884899431119, EPSILON); + EXPECT_NEAR(bonus[6].quat[2], -0.6379094464220707, EPSILON); + EXPECT_NEAR(bonus[6].quat[3], 0.1065166771202199, EPSILON); + EXPECT_NEAR(bonus[7].quat[0], -0.10483691088405202, EPSILON); + EXPECT_NEAR(bonus[7].quat[1], 0.5508895999584645, EPSILON); + EXPECT_NEAR(bonus[7].quat[2], -0.8250090480220789, EPSILON); + EXPECT_NEAR(bonus[7].quat[3], 0.06992811634525403, EPSILON); + EXPECT_NEAR(bonus[8].quat[0], 0.07777239911646, EPSILON); + EXPECT_NEAR(bonus[8].quat[1], -0.3724087549185288, EPSILON); + EXPECT_NEAR(bonus[8].quat[2], 0.9103052384821374, EPSILON); + EXPECT_NEAR(bonus[8].quat[3], -0.1631181963720798, EPSILON); + EXPECT_NEAR(bonus[9].quat[0], 0.16279109707978262, EPSILON); + EXPECT_NEAR(bonus[9].quat[1], 0.027148630125149613, EPSILON); + EXPECT_NEAR(bonus[9].quat[2], 0.9849325709665359, EPSILON); + EXPECT_NEAR(bonus[9].quat[3], -0.0516705065113425, EPSILON); + + auto num_bond = lmp->atom->num_bond; + auto bond_type = lmp->atom->bond_type; + auto bond_atom = lmp->atom->bond_atom; + auto id5p = lmp->atom->id5p; + + ASSERT_EQ(num_bond[GETIDX(1)], 1); + ASSERT_EQ(num_bond[GETIDX(2)], 1); + ASSERT_EQ(num_bond[GETIDX(3)], 1); + ASSERT_EQ(num_bond[GETIDX(4)], 1); + ASSERT_EQ(num_bond[GETIDX(5)], 0); + ASSERT_EQ(num_bond[GETIDX(6)], 1); + ASSERT_EQ(num_bond[GETIDX(7)], 1); + ASSERT_EQ(num_bond[GETIDX(8)], 1); + ASSERT_EQ(num_bond[GETIDX(9)], 1); + ASSERT_EQ(num_bond[GETIDX(10)], 0); + + ASSERT_EQ(bond_type[GETIDX(1)][0], 1); + ASSERT_EQ(bond_type[GETIDX(2)][0], 1); + ASSERT_EQ(bond_type[GETIDX(3)][0], 1); + ASSERT_EQ(bond_type[GETIDX(4)][0], 1); + ASSERT_EQ(bond_type[GETIDX(5)][0], 0); + ASSERT_EQ(bond_type[GETIDX(6)][0], 1); + ASSERT_EQ(bond_type[GETIDX(7)][0], 1); + ASSERT_EQ(bond_type[GETIDX(8)][0], 1); + ASSERT_EQ(bond_type[GETIDX(9)][0], 1); + ASSERT_EQ(bond_type[GETIDX(10)][0], 0); + + ASSERT_EQ(bond_atom[GETIDX(1)][0], 2); + ASSERT_EQ(bond_atom[GETIDX(2)][0], 3); + ASSERT_EQ(bond_atom[GETIDX(3)][0], 4); + ASSERT_EQ(bond_atom[GETIDX(4)][0], 5); + ASSERT_EQ(bond_atom[GETIDX(5)][0], 0); + ASSERT_EQ(bond_atom[GETIDX(6)][0], 7); + ASSERT_EQ(bond_atom[GETIDX(7)][0], 8); + ASSERT_EQ(bond_atom[GETIDX(8)][0], 9); + ASSERT_EQ(bond_atom[GETIDX(9)][0], 10); + ASSERT_EQ(bond_atom[GETIDX(10)][0], 0); + + ASSERT_EQ(id5p[GETIDX(1)], 2); + ASSERT_EQ(id5p[GETIDX(2)], 3); + ASSERT_EQ(id5p[GETIDX(3)], 4); + ASSERT_EQ(id5p[GETIDX(4)], 5); + ASSERT_EQ(id5p[GETIDX(5)], -1); + ASSERT_EQ(id5p[GETIDX(6)], 7); + ASSERT_EQ(id5p[GETIDX(7)], 8); + ASSERT_EQ(id5p[GETIDX(8)], 9); + ASSERT_EQ(id5p[GETIDX(9)], 10); + ASSERT_EQ(id5p[GETIDX(10)], -1); + + END_HIDE_OUTPUT(); + +} + } // namespace LAMMPS_NS int main(int argc, char **argv) From d3c45d33898b60d98aa4f9bc84a2aa85889a9832 Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Fri, 28 Jan 2022 17:58:00 +0000 Subject: [PATCH 11/13] Removed whitespace --- src/CG-DNA/pair_oxdna2_coaxstk.h | 9 +-------- src/CG-DNA/pair_oxdna2_dh.h | 4 +--- src/CG-DNA/pair_oxdna_coaxstk.h | 6 ------ src/CG-DNA/pair_oxdna_excv.h | 1 - src/CG-DNA/pair_oxdna_hbond.h | 8 -------- src/CG-DNA/pair_oxdna_stk.h | 3 --- src/CG-DNA/pair_oxdna_xstk.h | 8 -------- src/CG-DNA/pair_oxrna2_stk.h | 5 +---- src/CG-DNA/pair_oxrna2_xstk.h | 9 +-------- 9 files changed, 4 insertions(+), 49 deletions(-) diff --git a/src/CG-DNA/pair_oxdna2_coaxstk.h b/src/CG-DNA/pair_oxdna2_coaxstk.h index 7b376fa368..eecd369c38 100644 --- a/src/CG-DNA/pair_oxdna2_coaxstk.h +++ b/src/CG-DNA/pair_oxdna2_coaxstk.h @@ -46,23 +46,16 @@ class PairOxdna2Coaxstk : public Pair { double **k_cxst, **cut_cxst_0, **cut_cxst_c, **cut_cxst_lo, **cut_cxst_hi; double **cut_cxst_lc, **cut_cxst_hc, **b_cxst_lo, **b_cxst_hi; double **cutsq_cxst_hc; - double **a_cxst1, **theta_cxst1_0, **dtheta_cxst1_ast; double **b_cxst1, **dtheta_cxst1_c; - double **a_cxst4, **theta_cxst4_0, **dtheta_cxst4_ast; double **b_cxst4, **dtheta_cxst4_c; - double **a_cxst5, **theta_cxst5_0, **dtheta_cxst5_ast; double **b_cxst5, **dtheta_cxst5_c; - double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast; double **b_cxst6, **dtheta_cxst6_c; - double **AA_cxst1, **BB_cxst1; - - // per-atom arrays for local unit vectors - double **nx_xtrct, **nz_xtrct; + double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxdna2_dh.h b/src/CG-DNA/pair_oxdna2_dh.h index d7075f2988..dd299a197c 100644 --- a/src/CG-DNA/pair_oxdna2_dh.h +++ b/src/CG-DNA/pair_oxdna2_dh.h @@ -45,9 +45,7 @@ class PairOxdna2Dh : public Pair { protected: double **qeff_dh_pf, **kappa_dh; double **b_dh, **cut_dh_ast, **cutsq_dh_ast, **cut_dh_c, **cutsq_dh_c; - - // per-atom arrays for local unit vectors - double **nx_xtrct, **ny_xtrct, **nz_xtrct; + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxdna_coaxstk.h b/src/CG-DNA/pair_oxdna_coaxstk.h index 28813a034f..1868e9897e 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.h +++ b/src/CG-DNA/pair_oxdna_coaxstk.h @@ -47,22 +47,16 @@ class PairOxdnaCoaxstk : public Pair { double **k_cxst, **cut_cxst_0, **cut_cxst_c, **cut_cxst_lo, **cut_cxst_hi; double **cut_cxst_lc, **cut_cxst_hc, **b_cxst_lo, **b_cxst_hi; double **cutsq_cxst_hc; - double **a_cxst1, **theta_cxst1_0, **dtheta_cxst1_ast; double **b_cxst1, **dtheta_cxst1_c; - double **a_cxst4, **theta_cxst4_0, **dtheta_cxst4_ast; double **b_cxst4, **dtheta_cxst4_c; - double **a_cxst5, **theta_cxst5_0, **dtheta_cxst5_ast; double **b_cxst5, **dtheta_cxst5_c; - double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast; double **b_cxst6, **dtheta_cxst6_c; - double **a_cxst3p, **cosphi_cxst3p_ast, **b_cxst3p, **cosphi_cxst3p_c; double **a_cxst4p, **cosphi_cxst4p_ast, **b_cxst4p, **cosphi_cxst4p_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors virtual void allocate(); diff --git a/src/CG-DNA/pair_oxdna_excv.h b/src/CG-DNA/pair_oxdna_excv.h index 9abd97fb3d..c5bdc1a85a 100644 --- a/src/CG-DNA/pair_oxdna_excv.h +++ b/src/CG-DNA/pair_oxdna_excv.h @@ -55,7 +55,6 @@ class PairOxdnaExcv : public Pair { double **lj1_sb, **lj2_sb, **b_sb, **cut_sb_c, **cutsq_sb_c; double **epsilon_bb, **sigma_bb, **cut_bb_ast, **cutsq_bb_ast; double **lj1_bb, **lj2_bb, **b_bb, **cut_bb_c, **cutsq_bb_c; - double **nx, **ny, **nz; // per-atom arrays for local unit vectors virtual void allocate(); diff --git a/src/CG-DNA/pair_oxdna_hbond.h b/src/CG-DNA/pair_oxdna_hbond.h index 4afd9b8ad6..14bbfc5573 100644 --- a/src/CG-DNA/pair_oxdna_hbond.h +++ b/src/CG-DNA/pair_oxdna_hbond.h @@ -48,27 +48,19 @@ class PairOxdnaHbond : public Pair { double **epsilon_hb, **a_hb, **cut_hb_0, **cut_hb_c, **cut_hb_lo, **cut_hb_hi; double **cut_hb_lc, **cut_hb_hc, **b_hb_lo, **b_hb_hi, **shift_hb; double **cutsq_hb_hc; - double **a_hb1, **theta_hb1_0, **dtheta_hb1_ast; double **b_hb1, **dtheta_hb1_c; - double **a_hb2, **theta_hb2_0, **dtheta_hb2_ast; double **b_hb2, **dtheta_hb2_c; - double **a_hb3, **theta_hb3_0, **dtheta_hb3_ast; double **b_hb3, **dtheta_hb3_c; - double **a_hb4, **theta_hb4_0, **dtheta_hb4_ast; double **b_hb4, **dtheta_hb4_c; - double **a_hb7, **theta_hb7_0, **dtheta_hb7_ast; double **b_hb7, **dtheta_hb7_c; - double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast; double **b_hb8, **dtheta_hb8_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors - int seqdepflag; virtual void allocate(); diff --git a/src/CG-DNA/pair_oxdna_stk.h b/src/CG-DNA/pair_oxdna_stk.h index 5dbe9ea1e9..da2ad81673 100644 --- a/src/CG-DNA/pair_oxdna_stk.h +++ b/src/CG-DNA/pair_oxdna_stk.h @@ -59,10 +59,7 @@ class PairOxdnaStk : public Pair { double **b_st6, **dtheta_st6_c; double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c; double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors - - int seqdepflag; virtual void allocate(); diff --git a/src/CG-DNA/pair_oxdna_xstk.h b/src/CG-DNA/pair_oxdna_xstk.h index 7695effc89..04f9462fd4 100644 --- a/src/CG-DNA/pair_oxdna_xstk.h +++ b/src/CG-DNA/pair_oxdna_xstk.h @@ -47,28 +47,20 @@ class PairOxdnaXstk : public Pair { double **k_xst, **cut_xst_0, **cut_xst_c, **cut_xst_lo, **cut_xst_hi; double **cut_xst_lc, **cut_xst_hc, **b_xst_lo, **b_xst_hi; double **cutsq_xst_hc; - double **a_xst1, **theta_xst1_0, **dtheta_xst1_ast; double **b_xst1, **dtheta_xst1_c; - double **a_xst2, **theta_xst2_0, **dtheta_xst2_ast; double **b_xst2, **dtheta_xst2_c; - double **a_xst3, **theta_xst3_0, **dtheta_xst3_ast; double **b_xst3, **dtheta_xst3_c; - double **a_xst4, **theta_xst4_0, **dtheta_xst4_ast; double **b_xst4, **dtheta_xst4_c; - double **a_xst7, **theta_xst7_0, **dtheta_xst7_ast; double **b_xst7, **dtheta_xst7_c; - double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast; double **b_xst8, **dtheta_xst8_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors - virtual void allocate(); }; diff --git a/src/CG-DNA/pair_oxrna2_stk.h b/src/CG-DNA/pair_oxrna2_stk.h index e15614f126..0e7647894d 100644 --- a/src/CG-DNA/pair_oxrna2_stk.h +++ b/src/CG-DNA/pair_oxrna2_stk.h @@ -60,10 +60,7 @@ class PairOxrna2Stk : public Pair { double **b_st10, **dtheta_st10_c; double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c; double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c; - - // per-atom arrays for local unit vectors - double **nx_xtrct, **ny_xtrct, **nz_xtrct; - + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors int seqdepflag; virtual void allocate(); diff --git a/src/CG-DNA/pair_oxrna2_xstk.h b/src/CG-DNA/pair_oxrna2_xstk.h index 6ece8e3f32..8516b1c5d4 100644 --- a/src/CG-DNA/pair_oxrna2_xstk.h +++ b/src/CG-DNA/pair_oxrna2_xstk.h @@ -46,24 +46,17 @@ class PairOxrna2Xstk : public Pair { double **k_xst, **cut_xst_0, **cut_xst_c, **cut_xst_lo, **cut_xst_hi; double **cut_xst_lc, **cut_xst_hc, **b_xst_lo, **b_xst_hi; double **cutsq_xst_hc; - double **a_xst1, **theta_xst1_0, **dtheta_xst1_ast; double **b_xst1, **dtheta_xst1_c; - double **a_xst2, **theta_xst2_0, **dtheta_xst2_ast; double **b_xst2, **dtheta_xst2_c; - double **a_xst3, **theta_xst3_0, **dtheta_xst3_ast; double **b_xst3, **dtheta_xst3_c; - double **a_xst7, **theta_xst7_0, **dtheta_xst7_ast; double **b_xst7, **dtheta_xst7_c; - double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast; double **b_xst8, **dtheta_xst8_c; - - // per-atom arrays for local unit vectors - double **nx_xtrct, **nz_xtrct; + double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors virtual void allocate(); }; From 1d51f2d15177d87a0caa527e4e16eedd2fcb5ad6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 28 Jan 2022 16:57:58 -0500 Subject: [PATCH 12/13] must not access invalid/uninitialized data --- unittest/formats/test_atom_styles.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/unittest/formats/test_atom_styles.cpp b/unittest/formats/test_atom_styles.cpp index 5076e92ca3..7329c6d7a0 100644 --- a/unittest/formats/test_atom_styles.cpp +++ b/unittest/formats/test_atom_styles.cpp @@ -5220,23 +5220,19 @@ TEST_F(AtomStyleTest, oxdna) ASSERT_EQ(bond_type[GETIDX(2)][0], 1); ASSERT_EQ(bond_type[GETIDX(3)][0], 1); ASSERT_EQ(bond_type[GETIDX(4)][0], 1); - ASSERT_EQ(bond_type[GETIDX(5)][0], 0); ASSERT_EQ(bond_type[GETIDX(6)][0], 1); ASSERT_EQ(bond_type[GETIDX(7)][0], 1); ASSERT_EQ(bond_type[GETIDX(8)][0], 1); ASSERT_EQ(bond_type[GETIDX(9)][0], 1); - ASSERT_EQ(bond_type[GETIDX(10)][0], 0); ASSERT_EQ(bond_atom[GETIDX(1)][0], 2); ASSERT_EQ(bond_atom[GETIDX(2)][0], 3); ASSERT_EQ(bond_atom[GETIDX(3)][0], 4); ASSERT_EQ(bond_atom[GETIDX(4)][0], 5); - ASSERT_EQ(bond_atom[GETIDX(5)][0], 0); ASSERT_EQ(bond_atom[GETIDX(6)][0], 7); ASSERT_EQ(bond_atom[GETIDX(7)][0], 8); ASSERT_EQ(bond_atom[GETIDX(8)][0], 9); ASSERT_EQ(bond_atom[GETIDX(9)][0], 10); - ASSERT_EQ(bond_atom[GETIDX(10)][0], 0); ASSERT_EQ(id5p[GETIDX(1)], 2); ASSERT_EQ(id5p[GETIDX(2)], 3); From e8ce01079d3ad2aab45f635c6e9cd046278787f4 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 28 Jan 2022 17:20:32 -0500 Subject: [PATCH 13/13] whitespace --- src/CG-DNA/pair_oxdna_coaxstk.cpp | 4 ++-- src/CG-DNA/pair_oxrna2_stk.h | 2 +- src/CG-DNA/pair_oxrna2_xstk.h | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CG-DNA/pair_oxdna_coaxstk.cpp b/src/CG-DNA/pair_oxdna_coaxstk.cpp index e6cf59f072..49b25bb81b 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna_coaxstk.cpp @@ -158,7 +158,7 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) alist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - + // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv int dim; nx_xtrct = (double **) force->pair->extract("nx",dim); @@ -252,7 +252,7 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) // early rejection criterium if (f4t1) { - + az[0] = nz_xtrct[a][0]; az[1] = nz_xtrct[a][1]; az[2] = nz_xtrct[a][2]; diff --git a/src/CG-DNA/pair_oxrna2_stk.h b/src/CG-DNA/pair_oxrna2_stk.h index dc365f1c44..8530544fdc 100644 --- a/src/CG-DNA/pair_oxrna2_stk.h +++ b/src/CG-DNA/pair_oxrna2_stk.h @@ -60,7 +60,7 @@ class PairOxrna2Stk : public Pair { double **b_st10, **dtheta_st10_c; double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c; double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c; - double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors + double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors int seqdepflag; virtual void allocate(); diff --git a/src/CG-DNA/pair_oxrna2_xstk.h b/src/CG-DNA/pair_oxrna2_xstk.h index 90d7960f84..c73839b4a7 100644 --- a/src/CG-DNA/pair_oxrna2_xstk.h +++ b/src/CG-DNA/pair_oxrna2_xstk.h @@ -56,7 +56,7 @@ class PairOxrna2Xstk : public Pair { double **b_xst7, **dtheta_xst7_c; double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast; double **b_xst8, **dtheta_xst8_c; - double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors + double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors virtual void allocate(); };