diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 811164ff3..31e14adcf 100644 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -1,186 +1,200 @@ # Install/unInstall package files in LAMMPS # mode = 0/1/2 for uninstall/install/update mode=$1 # arg1 = file, arg2 = file it depends on action () { if (test $mode = 0) then rm -f ../$1 elif (! cmp -s $1 ../$1) then if (test -z "$2" || test -e ../$2) then cp $1 .. if (test $mode = 2) then echo " updating src/$1" fi fi elif (test -n "$2") then if (test ! -e ../$2) then rm -f ../$1 fi fi } # force rebuild of files with LMP_KOKKOS switch touch ../accelerator_kokkos.h touch ../memory.h # list of files with optional dependcies action angle_charmm_kokkos.cpp angle_charmm.cpp action angle_charmm_kokkos.h angle_charmm.h action angle_harmonic_kokkos.cpp angle_harmonic.cpp action angle_harmonic_kokkos.h angle_harmonic.h action atom_kokkos.cpp action atom_kokkos.h action atom_vec_angle_kokkos.cpp atom_vec_angle.cpp action atom_vec_angle_kokkos.h atom_vec_angle.h action atom_vec_atomic_kokkos.cpp action atom_vec_atomic_kokkos.h action atom_vec_bond_kokkos.cpp atom_vec_bond.cpp action atom_vec_bond_kokkos.h atom_vec_bond.h action atom_vec_charge_kokkos.cpp action atom_vec_charge_kokkos.h action atom_vec_full_kokkos.cpp atom_vec_full.cpp action atom_vec_full_kokkos.h atom_vec_full.h action atom_vec_kokkos.cpp action atom_vec_kokkos.h action atom_vec_molecular_kokkos.cpp atom_vec_molecular.cpp action atom_vec_molecular_kokkos.h atom_vec_molecular.h action bond_fene_kokkos.cpp bond_fene.cpp action bond_fene_kokkos.h bond_fene.h action bond_harmonic_kokkos.cpp bond_harmonic.cpp action bond_harmonic_kokkos.h bond_harmonic.h action comm_kokkos.cpp action comm_kokkos.h +action compute_temp_kokkos.cpp +action compute_temp_kokkos.h action dihedral_charmm_kokkos.cpp dihedral_charmm.cpp action dihedral_charmm_kokkos.h dihedral_charmm.h action dihedral_opls_kokkos.cpp dihedral_opls.cpp action dihedral_opls_kokkos.h dihedral_opls.h action domain_kokkos.cpp action domain_kokkos.h +action fix_deform_kokkos.cpp +action fix_deform_kokkos.h action fix_langevin_kokkos.cpp action fix_langevin_kokkos.h +action fix_nh_kokkos.cpp +action fix_nh_kokkos.h +action fix_nph_kokkos.cpp +action fix_nph_kokkos.h +action fix_npt_kokkos.cpp +action fix_npt_kokkos.h action fix_nve_kokkos.cpp action fix_nve_kokkos.h +action fix_nvt_kokkos.cpp +action fix_nvt_kokkos.h +action fix_wall_reflect_kokkos.cpp +action fix_wall_reflect_kokkos.h action improper_harmonic_kokkos.cpp improper_harmonic.cpp action improper_harmonic_kokkos.h improper_harmonic.h action kokkos.cpp action kokkos.h action kokkos_type.h action memory_kokkos.h action modify_kokkos.cpp action modify_kokkos.h action neigh_bond_kokkos.cpp action neigh_bond_kokkos.h action neigh_full_kokkos.h action neigh_list_kokkos.cpp action neigh_list_kokkos.h action neighbor_kokkos.cpp action neighbor_kokkos.h action pair_buck_coul_cut_kokkos.cpp action pair_buck_coul_cut_kokkos.h action pair_buck_coul_long_kokkos.cpp pair_buck_coul_long.cpp action pair_buck_coul_long_kokkos.h pair_buck_coul_long.h action pair_buck_kokkos.cpp action pair_buck_kokkos.h action pair_coul_cut_kokkos.cpp action pair_coul_cut_kokkos.h action pair_coul_debye_kokkos.cpp action pair_coul_debye_kokkos.h action pair_coul_dsf_kokkos.cpp action pair_coul_dsf_kokkos.h action pair_coul_long_kokkos.cpp pair_coul_long.cpp action pair_coul_long_kokkos.h pair_coul_long.h action pair_coul_wolf_kokkos.cpp action pair_coul_wolf_kokkos.h action pair_eam_kokkos.cpp pair_eam.cpp action pair_eam_kokkos.h pair_eam.h action pair_eam_alloy_kokkos.cpp pair_eam_alloy.cpp action pair_eam_alloy_kokkos.h pair_eam_alloy.h action pair_eam_fs_kokkos.cpp pair_eam_fs.cpp action pair_eam_fs_kokkos.h pair_eam_fs.h action pair_kokkos.h action pair_lj_charmm_coul_charmm_implicit_kokkos.cpp pair_lj_charmm_coul_charmm_implicit.cpp action pair_lj_charmm_coul_charmm_implicit_kokkos.h pair_lj_charmm_coul_charmm_implicit.h action pair_lj_charmm_coul_charmm_kokkos.cpp pair_lj_charmm_coul_charmm.cpp action pair_lj_charmm_coul_charmm_kokkos.h pair_lj_charmm_coul_charmm.h action pair_lj_charmm_coul_long_kokkos.cpp pair_lj_charmm_coul_long.cpp action pair_lj_charmm_coul_long_kokkos.h pair_lj_charmm_coul_long.h action pair_lj_class2_coul_cut_kokkos.cpp pair_lj_class2_coul_cut.cpp action pair_lj_class2_coul_cut_kokkos.h pair_lj_class2_coul_cut.h action pair_lj_class2_coul_long_kokkos.cpp pair_lj_class2_coul_long.cpp action pair_lj_class2_coul_long_kokkos.h pair_lj_class2_coul_long.h action pair_lj_class2_kokkos.cpp pair_lj_class2.cpp action pair_lj_class2_kokkos.h pair_lj_class2.h action pair_lj_cut_coul_cut_kokkos.cpp action pair_lj_cut_coul_cut_kokkos.h action pair_lj_cut_coul_debye_kokkos.cpp action pair_lj_cut_coul_debye_kokkos.h action pair_lj_cut_coul_dsf_kokkos.cpp action pair_lj_cut_coul_dsf_kokkos.h action pair_lj_cut_coul_long_kokkos.cpp pair_lj_cut_coul_long.cpp action pair_lj_cut_coul_long_kokkos.h pair_lj_cut_coul_long.h action pair_lj_cut_kokkos.cpp action pair_lj_cut_kokkos.h action pair_lj_expand_kokkos.cpp action pair_lj_expand_kokkos.h action pair_lj_gromacs_coul_gromacs_kokkos.cpp action pair_lj_gromacs_coul_gromacs_kokkos.h action pair_lj_gromacs_kokkos.cpp action pair_lj_gromacs_kokkos.h action pair_lj_sdk_kokkos.cpp pair_lj_sdk.cpp action pair_lj_sdk_kokkos.h pair_lj_sdk.h action pair_sw_kokkos.cpp pair_sw.cpp action pair_sw_kokkos.h pair_sw.h action pair_table_kokkos.cpp action pair_table_kokkos.h action pair_tersoff_kokkos.cpp pair_tersoff.cpp action pair_tersoff_kokkos.h pair_tersoff.h action pair_tersoff_mod_kokkos.cpp pair_tersoff_mod.cpp action pair_tersoff_mod_kokkos.h pair_tersoff_mod.h action pair_tersoff_zbl_kokkos.cpp pair_tersoff_zbl.cpp action pair_tersoff_zbl_kokkos.h pair_tersoff_zbl.h action verlet_kokkos.cpp action verlet_kokkos.h # edit 2 Makefile.package files to include/exclude package info if (test $1 = 1) then if (test -e ../Makefile.package) then sed -i -e 's/[^ \t]*kokkos[^ \t]* //g' ../Makefile.package sed -i -e 's/[^ \t]*KOKKOS[^ \t]* //g' ../Makefile.package sed -i -e 's|^PKG_INC =[ \t]*|&-DLMP_KOKKOS |' ../Makefile.package # sed -i -e 's|^PKG_PATH =[ \t]*|&-L..\/..\/lib\/kokkos\/core\/src |' ../Makefile.package sed -i -e 's|^PKG_CPP_DEPENDS =[ \t]*|&$(KOKKOS_CPP_DEPENDS) |' ../Makefile.package sed -i -e 's|^PKG_LIB =[ \t]*|&$(KOKKOS_LIBS) |' ../Makefile.package sed -i -e 's|^PKG_LINK_DEPENDS =[ \t]*|&$(KOKKOS_LINK_DEPENDS) |' ../Makefile.package sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(KOKKOS_CPPFLAGS) $(KOKKOS_CXXFLAGS) |' ../Makefile.package sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(KOKKOS_LDFLAGS) |' ../Makefile.package # sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(kokkos_SYSPATH) |' ../Makefile.package fi if (test -e ../Makefile.package.settings) then sed -i -e '/CXX\ =\ \$(CC)/d' ../Makefile.package.settings sed -i -e '/^include.*kokkos.*$/d' ../Makefile.package.settings # multiline form needed for BSD sed on Macs sed -i -e '4 i \CXX = $(CC)' ../Makefile.package.settings sed -i -e '5 i \include ..\/..\/lib\/kokkos\/Makefile.kokkos' ../Makefile.package.settings fi elif (test $1 = 0) then if (test -e ../Makefile.package) then sed -i -e 's/[^ \t]*kokkos[^ \t]* //g' ../Makefile.package sed -i -e 's/[^ \t]*KOKKOS[^ \t]* //g' ../Makefile.package fi if (test -e ../Makefile.package.settings) then sed -i -e '/CXX\ =\ \$(CC)/d' ../Makefile.package.settings sed -i -e '/^include.*kokkos.*$/d' ../Makefile.package.settings fi fi diff --git a/src/KOKKOS/domain_kokkos.cpp b/src/KOKKOS/domain_kokkos.cpp index 11c2147a0..5bb796b0f 100644 --- a/src/KOKKOS/domain_kokkos.cpp +++ b/src/KOKKOS/domain_kokkos.cpp @@ -1,207 +1,418 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "domain_kokkos.h" #include "atom_kokkos.h" #include "atom_masks.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ DomainKokkos::DomainKokkos(LAMMPS *lmp) : Domain(lmp) {} /* ---------------------------------------------------------------------- */ void DomainKokkos::init() { atomKK = (AtomKokkos *) atom; Domain::init(); } /* ---------------------------------------------------------------------- */ template struct DomainPBCFunctor { typedef DeviceType device_type; double lo[3],hi[3],period[3]; typename ArrayTypes::t_x_array x; typename ArrayTypes::t_v_array v; typename ArrayTypes::t_int_1d mask; typename ArrayTypes::t_imageint_1d image; int deform_groupbit; double h_rate[6]; int xperiodic,yperiodic,zperiodic; DomainPBCFunctor(double* _lo, double* _hi, double* _period, DAT::tdual_x_array _x, DAT::tdual_v_array _v, DAT::tdual_int_1d _mask, DAT::tdual_imageint_1d _image, int _deform_groupbit, double* _h_rate, int _xperiodic, int _yperiodic, int _zperiodic): x(_x.view()), v(_v.view()), mask(_mask.view()), image(_image.view()), deform_groupbit(_deform_groupbit), xperiodic(_xperiodic), yperiodic(_yperiodic), zperiodic(_zperiodic){ lo[0]=_lo[0]; lo[1]=_lo[1]; lo[2]=_lo[2]; hi[0]=_hi[0]; hi[1]=_hi[1]; hi[2]=_hi[2]; period[0]=_period[0]; period[1]=_period[1]; period[2]=_period[2]; h_rate[0]=_h_rate[0]; h_rate[1]=_h_rate[1]; h_rate[2]=_h_rate[2]; h_rate[3]=_h_rate[3]; h_rate[4]=_h_rate[4]; h_rate[5]=_h_rate[5]; } KOKKOS_INLINE_FUNCTION void operator() (const int &i) const { if (PERIODIC && xperiodic) { if (x(i,0) < lo[0]) { x(i,0) += period[0]; if (DEFORM_VREMAP && (mask[i] & deform_groupbit)) v(i,0) += h_rate[0]; imageint idim = image[i] & IMGMASK; const int otherdims = image[i] ^ idim; idim--; idim &= IMGMASK; image[i] = otherdims | idim; } if (x(i,0) >= hi[0]) { x(i,0) -= period[0]; x(i,0) = MAX(x(i,0),lo[0]); if (DEFORM_VREMAP && (mask[i] & deform_groupbit)) v(i,0) -= h_rate[0]; imageint idim = image[i] & IMGMASK; const int otherdims = image[i] ^ idim; idim++; idim &= IMGMASK; image[i] = otherdims | idim; } } if (PERIODIC && yperiodic) { if (x(i,1) < lo[1]) { x(i,1) += period[1]; if (DEFORM_VREMAP && (mask[i] & deform_groupbit)) { v(i,0) += h_rate[5]; v(i,1) += h_rate[1]; } imageint idim = (image[i] >> IMGBITS) & IMGMASK; const imageint otherdims = image[i] ^ (idim << IMGBITS); idim--; idim &= IMGMASK; image[i] = otherdims | (idim << IMGBITS); } if (x(i,1) >= hi[1]) { x(i,1) -= period[1]; x(i,1) = MAX(x(i,1),lo[1]); if (DEFORM_VREMAP && (mask[i] & deform_groupbit)) { v(i,0) -= h_rate[5]; v(i,1) -= h_rate[1]; } imageint idim = (image[i] >> IMGBITS) & IMGMASK; const imageint otherdims = image[i] ^ (idim << IMGBITS); idim++; idim &= IMGMASK; image[i] = otherdims | (idim << IMGBITS); } } if (PERIODIC && zperiodic) { if (x(i,2) < lo[2]) { x(i,2) += period[2]; if (DEFORM_VREMAP && (mask[i] & deform_groupbit)) { v(i,0) += h_rate[4]; v(i,1) += h_rate[3]; v(i,2) += h_rate[2]; } imageint idim = image[i] >> IMG2BITS; const imageint otherdims = image[i] ^ (idim << IMG2BITS); idim--; idim &= IMGMASK; image[i] = otherdims | (idim << IMG2BITS); } if (x(i,2) >= hi[2]) { x(i,2) -= period[2]; x(i,2) = MAX(x(i,2),lo[2]); if (DEFORM_VREMAP && (mask[i] & deform_groupbit)) { v(i,0) -= h_rate[4]; v(i,1) -= h_rate[3]; v(i,2) -= h_rate[2]; } imageint idim = image[i] >> IMG2BITS; const imageint otherdims = image[i] ^ (idim << IMG2BITS); idim++; idim &= IMGMASK; image[i] = otherdims | (idim << IMG2BITS); } } } }; /* ---------------------------------------------------------------------- enforce PBC and modify box image flags for each atom called every reneighboring and by other commands that change atoms resulting coord must satisfy lo <= coord < hi MAX is important since coord - prd < lo can happen when coord = hi if fix deform, remap velocity of fix group atoms by box edge velocities for triclinic, atoms must be in lamda coords (0-1) before pbc is called image = 10 bits for each dimension increment/decrement in wrap-around fashion ------------------------------------------------------------------------- */ void DomainKokkos::pbc() { double *lo,*hi,*period; int nlocal = atomKK->nlocal; if (triclinic == 0) { lo = boxlo; hi = boxhi; period = prd; } else { lo = boxlo_lamda; hi = boxhi_lamda; period = prd_lamda; } atomKK->sync(Device,X_MASK|V_MASK|MASK_MASK|IMAGE_MASK); atomKK->modified(Device,X_MASK|V_MASK); if (xperiodic || yperiodic || zperiodic) { if (deform_vremap) { DomainPBCFunctor f(lo,hi,period, atomKK->k_x,atomKK->k_v,atomKK->k_mask,atomKK->k_image, deform_groupbit,h_rate,xperiodic,yperiodic,zperiodic); Kokkos::parallel_for(nlocal,f); } else { DomainPBCFunctor f(lo,hi,period, atomKK->k_x,atomKK->k_v,atomKK->k_mask,atomKK->k_image, deform_groupbit,h_rate,xperiodic,yperiodic,zperiodic); Kokkos::parallel_for(nlocal,f); } } else { if (deform_vremap) { DomainPBCFunctor f(lo,hi,period, atomKK->k_x,atomKK->k_v,atomKK->k_mask,atomKK->k_image, deform_groupbit,h_rate,xperiodic,yperiodic,zperiodic); Kokkos::parallel_for(nlocal,f); } else { DomainPBCFunctor f(lo,hi,period, atomKK->k_x,atomKK->k_v,atomKK->k_mask,atomKK->k_image, deform_groupbit,h_rate,xperiodic,yperiodic,zperiodic); Kokkos::parallel_for(nlocal,f); } } LMPDeviceType::fence(); } +/* ---------------------------------------------------------------------- + remap all points into the periodic box no matter how far away + adjust 3 image flags encoded in image accordingly + resulting coord must satisfy lo <= coord < hi + MAX is important since coord - prd < lo can happen when coord = hi + for triclinic, point is converted to lamda coords (0-1) before doing remap + image = 10 bits for each dimension + increment/decrement in wrap-around fashion +------------------------------------------------------------------------- */ + +void DomainKokkos::remap_all() +{ + atomKK->sync(Device,X_MASK | IMAGE_MASK); + atomKK->modified(Device,X_MASK | IMAGE_MASK); + + x = atomKK->k_x.view(); + image = atomKK->k_image.view(); + int nlocal = atomKK->nlocal; + + if (triclinic == 0) { + for (int i=0; i<3; i++) { + lo[i] = boxlo[i]; + hi[i] = boxhi[i]; + period[i] = prd[i]; + } + } else { + for (int i=0; i<3; i++) { + lo[i] = boxlo_lamda[i]; + hi[i] = boxhi_lamda[i]; + period[i] = prd_lamda[i]; + } + x2lamda(nlocal); + } + + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); + LMPDeviceType::fence(); + copymode = 0; + + if (triclinic) lamda2x(nlocal); +} + +KOKKOS_INLINE_FUNCTION +void DomainKokkos::operator()(TagDomain_remap_all, const int &i) const { + imageint idim,otherdims; + if (xperiodic) { + while (x(i,0) < lo[0]) { + x(i,0) += period[0]; + idim = image[i] & IMGMASK; + otherdims = image[i] ^ idim; + idim--; + idim &= IMGMASK; + image[i] = otherdims | idim; + } + while (x(i,0) >= hi[0]) { + x(i,0) -= period[0]; + idim = image[i] & IMGMASK; + otherdims = image[i] ^ idim; + idim++; + idim &= IMGMASK; + image[i] = otherdims | idim; + } + x(i,0) = MAX(x(i,0),lo[0]); + } + + if (yperiodic) { + while (x(i,1) < lo[1]) { + x(i,1) += period[1]; + idim = (image[i] >> IMGBITS) & IMGMASK; + otherdims = image[i] ^ (idim << IMGBITS); + idim--; + idim &= IMGMASK; + image[i] = otherdims | (idim << IMGBITS); + } + while (x(i,1) >= hi[1]) { + x(i,1) -= period[1]; + idim = (image[i] >> IMGBITS) & IMGMASK; + otherdims = image[i] ^ (idim << IMGBITS); + idim++; + idim &= IMGMASK; + image[i] = otherdims | (idim << IMGBITS); + } + x(i,1) = MAX(x(i,1),lo[1]); + } + + if (zperiodic) { + while (x(i,2) < lo[2]) { + x(i,2) += period[2]; + idim = image[i] >> IMG2BITS; + otherdims = image[i] ^ (idim << IMG2BITS); + idim--; + idim &= IMGMASK; + image[i] = otherdims | (idim << IMG2BITS); + } + while (x(i,2) >= hi[2]) { + x(i,2) -= period[2]; + idim = image[i] >> IMG2BITS; + otherdims = image[i] ^ (idim << IMG2BITS); + idim++; + idim &= IMGMASK; + image[i] = otherdims | (idim << IMG2BITS); + } + x(i,2) = MAX(x(i,2),lo[2]); + } +} + +/* ---------------------------------------------------------------------- + adjust image flags due to triclinic box flip + flip operation is changing box vectors A,B,C to new A',B',C' + A' = A (A does not change) + B' = B + mA (B shifted by A) + C' = C + pB + nA (C shifted by B and/or A) + this requires the image flags change from (a,b,c) to (a',b',c') + so that x_unwrap for each atom is same before/after + x_unwrap_before = xlocal + aA + bB + cC + x_unwrap_after = xlocal + a'A' + b'B' + c'C' + this requires: + c' = c + b' = b - cp + a' = a - (b-cp)m - cn = a - b'm - cn + in other words, for xy flip, change in x flag depends on current y flag + this is b/c the xy flip dramatically changes which tiled image of + simulation box an unwrapped point maps to +------------------------------------------------------------------------- */ + +void DomainKokkos::image_flip(int m_in, int n_in, int p_in) +{ + m_flip = m_in; + n_flip = n_in; + p_flip = p_in; + + atomKK->sync(Device,IMAGE_MASK); + atomKK->modified(Device,IMAGE_MASK); + + image = atomKK->k_image.view(); + int nlocal = atomKK->nlocal; + + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); + LMPDeviceType::fence(); + copymode = 0; +} + +KOKKOS_INLINE_FUNCTION +void DomainKokkos::operator()(TagDomain_image_flip, const int &i) const { + int xbox = (image[i] & IMGMASK) - IMGMAX; + int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + int zbox = (image[i] >> IMG2BITS) - IMGMAX; + + ybox -= p_flip*zbox; + xbox -= m_flip*ybox + n_flip*zbox; + + image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) | + (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | + (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS); +} + +/* ---------------------------------------------------------------------- + convert triclinic 0-1 lamda coords to box coords for all N atoms + x = H lamda + x0; +------------------------------------------------------------------------- */ + +void DomainKokkos::lamda2x(int n) +{ + atomKK->sync(Device,X_MASK); + atomKK->modified(Device,X_MASK); + + x = atomKK->k_x.view(); + + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,n),*this); + LMPDeviceType::fence(); + copymode = 0; +} + +KOKKOS_INLINE_FUNCTION +void DomainKokkos::operator()(TagDomain_lamda2x, const int &i) const { + x(i,0) = h[0]*x(i,0) + h[5]*x(i,1) + h[4]*x(i,2) + boxlo[0]; + x(i,1) = h[1]*x(i,1) + h[3]*x(i,2) + boxlo[1]; + x(i,2) = h[2]*x(i,2) + boxlo[2]; +} + +/* ---------------------------------------------------------------------- + convert box coords to triclinic 0-1 lamda coords for all N atoms + lamda = H^-1 (x - x0) +------------------------------------------------------------------------- */ + +void DomainKokkos::x2lamda(int n) +{ + atomKK->sync(Device,X_MASK); + atomKK->modified(Device,X_MASK); + + x = atomKK->k_x.view(); + + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,n),*this); + LMPDeviceType::fence(); + copymode = 0; +} + +KOKKOS_INLINE_FUNCTION +void DomainKokkos::operator()(TagDomain_x2lamda, const int &i) const { + F_FLOAT delta[3]; + delta[0] = x(i,0) - boxlo[0]; + delta[1] = x(i,1) - boxlo[1]; + delta[2] = x(i,2) - boxlo[2]; + + x(i,0) = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2]; + x(i,1) = h_inv[1]*delta[1] + h_inv[3]*delta[2]; + x(i,2) = h_inv[2]*delta[2]; +} diff --git a/src/KOKKOS/domain_kokkos.h b/src/KOKKOS/domain_kokkos.h index 70ac7c974..c39b53a6c 100644 --- a/src/KOKKOS/domain_kokkos.h +++ b/src/KOKKOS/domain_kokkos.h @@ -1,38 +1,65 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifndef LMP_DOMAIN_KOKKOS_H #define LMP_DOMAIN_KOKKOS_H #include "domain.h" #include "kokkos_type.h" namespace LAMMPS_NS { +struct TagDomain_remap_all{}; +struct TagDomain_image_flip{}; +struct TagDomain_lamda2x{}; +struct TagDomain_x2lamda{}; + class DomainKokkos : public Domain { public: - - DomainKokkos(class LAMMPS *); ~DomainKokkos() {} void init(); void pbc(); + void remap_all(); + void image_flip(int, int, int); + void x2lamda(int); + void lamda2x(int); + + int closest_image(const int, int) const; + + KOKKOS_INLINE_FUNCTION + void operator()(TagDomain_remap_all, const int&) const; + + KOKKOS_INLINE_FUNCTION + void operator()(TagDomain_image_flip, const int&) const; + + KOKKOS_INLINE_FUNCTION + void operator()(TagDomain_lamda2x, const int&) const; + + KOKKOS_INLINE_FUNCTION + void operator()(TagDomain_x2lamda, const int&) const; + + private: + double lo[3],hi[3],period[3]; + int n_flip, m_flip, p_flip; + ArrayTypes::t_x_array x; + ArrayTypes::t_imageint_1d image; }; } #endif /* ERROR/WARNING messages: */ diff --git a/src/compute.cpp b/src/compute.cpp index 54dc57809..4c5da8bc9 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -1,233 +1,241 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "mpi.h" #include "stdlib.h" #include "string.h" #include "ctype.h" #include "compute.h" #include "atom.h" #include "domain.h" #include "comm.h" #include "group.h" #include "modify.h" #include "fix.h" #include "atom_masks.h" #include "memory.h" #include "error.h" #include "force.h" using namespace LAMMPS_NS; #define DELTA 4 #define BIG MAXTAGINT // allocate space for static class instance variable and initialize it int Compute::instance_total = 0; /* ---------------------------------------------------------------------- */ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) { instance_me = instance_total++; if (narg < 3) error->all(FLERR,"Illegal compute command"); // compute ID, group, and style // ID must be all alphanumeric chars or underscores int n = strlen(arg[0]) + 1; id = new char[n]; strcpy(id,arg[0]); for (int i = 0; i < n-1; i++) if (!isalnum(id[i]) && id[i] != '_') error->all(FLERR, "Compute ID must be alphanumeric or underscore characters"); igroup = group->find(arg[1]); if (igroup == -1) error->all(FLERR,"Could not find compute group ID"); groupbit = group->bitmask[igroup]; n = strlen(arg[2]) + 1; style = new char[n]; strcpy(style,arg[2]); // set child class defaults scalar_flag = vector_flag = array_flag = 0; peratom_flag = local_flag = 0; size_vector_variable = size_array_rows_variable = 0; tempflag = pressflag = peflag = 0; pressatomflag = peatomflag = 0; create_attribute = 0; tempbias = 0; timeflag = 0; comm_forward = comm_reverse = 0; dynamic = 0; dynamic_group_allow = 1; cudable = 0; invoked_scalar = invoked_vector = invoked_array = -1; invoked_peratom = invoked_local = -1; invoked_flag = 0; // set modify defaults extra_dof = domain->dimension; dynamic_user = 0; fix_dof = 0; // setup list of timesteps ntime = maxtime = 0; tlist = NULL; // data masks datamask = ALL_MASK; datamask_ext = ALL_MASK; + execution_space = Host; + datamask_read = ALL_MASK; + datamask_modify = ALL_MASK; + + copymode = 0; + // force init to zero in case these are used as logicals vector = vector_atom = vector_local = NULL; array = array_atom = array_local = NULL; } /* ---------------------------------------------------------------------- */ Compute::~Compute() { + if (copymode) return; + delete [] id; delete [] style; memory->destroy(tlist); } /* ---------------------------------------------------------------------- */ void Compute::modify_params(int narg, char **arg) { if (narg == 0) error->all(FLERR,"Illegal compute_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"extra") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command"); extra_dof = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"dynamic") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command"); if (strcmp(arg[iarg+1],"no") == 0) dynamic_user = 0; else if (strcmp(arg[iarg+1],"yes") == 0) dynamic_user = 1; else error->all(FLERR,"Illegal compute_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"thermo") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command"); if (strcmp(arg[iarg+1],"no") == 0) thermoflag = 0; else if (strcmp(arg[iarg+1],"yes") == 0) thermoflag = 1; else error->all(FLERR,"Illegal compute_modify command"); iarg += 2; } else error->all(FLERR,"Illegal compute_modify command"); } } /* ---------------------------------------------------------------------- calculate adjustment in DOF due to fixes ------------------------------------------------------------------------- */ void Compute::adjust_dof_fix() { Fix **fix = modify->fix; int nfix = modify->nfix; fix_dof = 0; for (int i = 0; i < nfix; i++) if (fix[i]->dof_flag) fix_dof += fix[i]->dof(igroup); } /* ---------------------------------------------------------------------- reset extra_dof to its default value ------------------------------------------------------------------------- */ void Compute::reset_extra_dof() { extra_dof = domain->dimension; } /* ---------------------------------------------------------------------- */ void Compute::reset_extra_compute_fix(const char *) { error->all(FLERR, "Compute does not allow an extra compute or fix to be reset"); } /* ---------------------------------------------------------------------- add ntimestep to list of timesteps the compute will be called on do not add if already in list search from top downward, since list of times is in decreasing order ------------------------------------------------------------------------- */ void Compute::addstep(bigint ntimestep) { // i = location in list to insert ntimestep int i; for (i = ntime-1; i >= 0; i--) { if (ntimestep == tlist[i]) return; if (ntimestep < tlist[i]) break; } i++; // extend list as needed if (ntime == maxtime) { maxtime += DELTA; memory->grow(tlist,maxtime,"compute:tlist"); } // move remainder of list upward and insert ntimestep for (int j = ntime-1; j >= i; j--) tlist[j+1] = tlist[j]; tlist[i] = ntimestep; ntime++; } /* ---------------------------------------------------------------------- return 1/0 if ntimestep is or is not in list of calling timesteps if value(s) on top of list are less than ntimestep, delete them search from top downward, since list of times is in decreasing order ------------------------------------------------------------------------- */ int Compute::matchstep(bigint ntimestep) { for (int i = ntime-1; i >= 0; i--) { if (ntimestep < tlist[i]) return 0; if (ntimestep == tlist[i]) return 1; if (ntimestep > tlist[i]) ntime--; } return 0; } /* ---------------------------------------------------------------------- clean out list of timesteps to call the compute on ------------------------------------------------------------------------- */ void Compute::clearstep() { ntime = 0; } diff --git a/src/compute.h b/src/compute.h index 732ced15d..29569a008 100644 --- a/src/compute.h +++ b/src/compute.h @@ -1,200 +1,207 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifndef LMP_COMPUTE_H #define LMP_COMPUTE_H #include "pointers.h" namespace LAMMPS_NS { class Compute : protected Pointers { public: static int instance_total; // # of Compute classes ever instantiated char *id,*style; int igroup,groupbit; double scalar; // computed global scalar double *vector; // computed global vector double **array; // computed global array double *vector_atom; // computed per-atom vector double **array_atom; // computed per-atom array double *vector_local; // computed local vector double **array_local; // computed local array int scalar_flag; // 0/1 if compute_scalar() function exists int vector_flag; // 0/1 if compute_vector() function exists int array_flag; // 0/1 if compute_array() function exists int size_vector; // length of global vector int size_array_rows; // rows in global array int size_array_cols; // columns in global array int size_vector_variable; // 1 if vec length is unknown in advance int size_array_rows_variable; // 1 if array rows is unknown in advance int peratom_flag; // 0/1 if compute_peratom() function exists int size_peratom_cols; // 0 = vector, N = columns in peratom array int local_flag; // 0/1 if compute_local() function exists int size_local_rows; // rows in local vector or array int size_local_cols; // 0 = vector, N = columns in local array int extscalar; // 0/1 if global scalar is intensive/extensive int extvector; // 0/1/-1 if global vector is all int/ext/extlist int *extlist; // list of 0/1 int/ext for each vec component int extarray; // 0/1 if global array is all intensive/extensive int tempflag; // 1 if Compute can be used as temperature // must have both compute_scalar, compute_vector int pressflag; // 1 if Compute can be used as pressure (uses virial) // must have both compute_scalar, compute_vector int pressatomflag; // 1 if Compute calculates per-atom virial int peflag; // 1 if Compute calculates PE (uses Force energies) int peatomflag; // 1 if Compute calculates per-atom PE int create_attribute; // 1 if compute stores attributes that need // setting when a new atom is created int tempbias; // 0/1 if Compute temp includes self/extra bias int timeflag; // 1 if Compute stores list of timesteps it's called on int ntime; // # of entries in time list int maxtime; // max # of entries time list can hold bigint *tlist; // list of timesteps the Compute is called on int invoked_flag; // non-zero if invoked or accessed this step, 0 if not bigint invoked_scalar; // last timestep on which compute_scalar() was invoked bigint invoked_vector; // ditto for compute_vector() bigint invoked_array; // ditto for compute_array() bigint invoked_peratom; // ditto for compute_peratom() bigint invoked_local; // ditto for compute_local() double dof; // degrees-of-freedom for temperature int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 unsigned int datamask; unsigned int datamask_ext; + // KOKKOS host/device flag and data masks + + ExecutionSpace execution_space; + unsigned int datamask_read,datamask_modify; + + int copymode; + int cudable; // 1 if compute is CUDA-enabled Compute(class LAMMPS *, int, char **); virtual ~Compute(); void modify_params(int, char **); void reset_extra_dof(); virtual void init() = 0; virtual void init_list(int, class NeighList *) {} virtual void setup() {} virtual double compute_scalar() {return 0.0;} virtual void compute_vector() {} virtual void compute_array() {} virtual void compute_peratom() {} virtual void compute_local() {} virtual void set_arrays(int) {} virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} virtual void dof_remove_pre() {} virtual int dof_remove(int) {return 0;} virtual void remove_bias(int, double *) {} virtual void remove_bias_all() {} virtual void reapply_bias_all() {} virtual void restore_bias(int, double *) {} virtual void restore_bias_all() {} virtual void reset_extra_compute_fix(const char *); virtual void lock_enable() {} virtual void lock_disable() {} virtual int lock_length() {return 0;} virtual void lock(class Fix *, bigint, bigint) {} virtual void unlock(class Fix *) {} void addstep(bigint); int matchstep(bigint); void clearstep(); virtual double memory_usage() {return 0.0;} virtual void pair_tally_callback(int, int, int, int, double, double, double, double, double, double) {} virtual int unsigned data_mask() {return datamask;} virtual int unsigned data_mask_ext() {return datamask_ext;} protected: int instance_me; // which Compute class instantiation I am double natoms_temp; // # of atoms used for temperature calculation int extra_dof; // extra DOF for temperature computes int fix_dof; // DOF due to fixes int dynamic; // recount atoms for temperature computes int dynamic_user; // user request for temp compute to be dynamic int thermoflag; // 1 if include fix PE for PE computes double vbias[3]; // stored velocity bias for one atom double **vbiasall; // stored velocity bias for all atoms int maxbias; // size of vbiasall array inline int sbmask(int j) { return j >> SBBITS & 3; } // union data struct for packing 32-bit and 64-bit ints into double bufs // see atom_vec.h for documentation union ubuf { double d; int64_t i; ubuf(double arg) : d(arg) {} ubuf(int64_t arg) : i(arg) {} ubuf(int arg) : i(arg) {} }; // private methods void adjust_dof_fix(); }; } #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Compute ID must be alphanumeric or underscore characters Self-explanatory. E: Could not find compute group ID Self-explanatory. E: Compute does not allow an extra compute or fix to be reset This is an internal LAMMPS error. Please report it to the developers. */ diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index 4a830d96a..6748bea4e 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -1,137 +1,138 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "mpi.h" #include "string.h" #include "compute_temp.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "comm.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTemp::ComputeTemp(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTemp::~ComputeTemp() { - delete [] vector; + if (!copymode) + delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTemp::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTemp::dof_compute() { adjust_dof_fix(); natoms_temp = group->count(igroup); dof = domain->dimension * natoms_temp; dof -= extra_dof + fix_dof; if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTemp::compute_scalar() { invoked_scalar = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); if (dof < 0.0 && natoms_temp > 0.0) error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTemp::compute_vector() { int i; invoked_vector = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; t[3] += massone * v[i][0]*v[i][1]; t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } diff --git a/src/compute_temp.h b/src/compute_temp.h index cfe5012f8..71175b677 100644 --- a/src/compute_temp.h +++ b/src/compute_temp.h @@ -1,55 +1,55 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef COMPUTE_CLASS ComputeStyle(temp,ComputeTemp) #else #ifndef LMP_COMPUTE_TEMP_H #define LMP_COMPUTE_TEMP_H #include "compute.h" namespace LAMMPS_NS { class ComputeTemp : public Compute { public: ComputeTemp(class LAMMPS *, int, char **); virtual ~ComputeTemp(); void init() {} void setup(); - double compute_scalar(); - void compute_vector(); + virtual double compute_scalar(); + virtual void compute_vector(); protected: double tfactor; virtual void dof_compute(); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. */ diff --git a/src/domain.cpp b/src/domain.cpp index 66d5f5936..b19a7c508 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -1,1868 +1,1872 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author (triclinic) : Pieter in 't Veld (SNL) ------------------------------------------------------------------------- */ #include "mpi.h" #include "stdlib.h" #include "string.h" #include "stdio.h" #include "math.h" #include "domain.h" #include "style_region.h" #include "atom.h" #include "atom_vec.h" #include "molecule.h" #include "force.h" #include "kspace.h" #include "update.h" #include "modify.h" #include "fix.h" #include "fix_deform.h" #include "region.h" #include "lattice.h" #include "comm.h" #include "output.h" #include "thermo.h" #include "universe.h" #include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp enum{IGNORE,WARN,ERROR}; // same as thermo.cpp enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files #define BIG 1.0e20 #define SMALL 1.0e-4 #define DELTAREGION 4 #define BONDSTRETCH 1.1 /* ---------------------------------------------------------------------- default is periodic ------------------------------------------------------------------------- */ Domain::Domain(LAMMPS *lmp) : Pointers(lmp) { box_exist = 0; dimension = 3; nonperiodic = 0; xperiodic = yperiodic = zperiodic = 1; periodicity[0] = xperiodic; periodicity[1] = yperiodic; periodicity[2] = zperiodic; boundary[0][0] = boundary[0][1] = 0; boundary[1][0] = boundary[1][1] = 0; boundary[2][0] = boundary[2][1] = 0; minxlo = minxhi = 0.0; minylo = minyhi = 0.0; minzlo = minzhi = 0.0; triclinic = 0; tiltsmall = 1; boxlo[0] = boxlo[1] = boxlo[2] = -0.5; boxhi[0] = boxhi[1] = boxhi[2] = 0.5; xy = xz = yz = 0.0; h[3] = h[4] = h[5] = 0.0; h_inv[3] = h_inv[4] = h_inv[5] = 0.0; h_rate[0] = h_rate[1] = h_rate[2] = h_rate[3] = h_rate[4] = h_rate[5] = 0.0; h_ratelo[0] = h_ratelo[1] = h_ratelo[2] = 0.0; prd_lamda[0] = prd_lamda[1] = prd_lamda[2] = 1.0; prd_half_lamda[0] = prd_half_lamda[1] = prd_half_lamda[2] = 0.5; boxlo_lamda[0] = boxlo_lamda[1] = boxlo_lamda[2] = 0.0; boxhi_lamda[0] = boxhi_lamda[1] = boxhi_lamda[2] = 1.0; lattice = NULL; char **args = new char*[2]; args[0] = (char *) "none"; args[1] = (char *) "1.0"; set_lattice(2,args); delete [] args; nregion = maxregion = 0; regions = NULL; + + copymode = 0; } /* ---------------------------------------------------------------------- */ Domain::~Domain() { + if (copymode) return; + delete lattice; for (int i = 0; i < nregion; i++) delete regions[i]; memory->sfree(regions); } /* ---------------------------------------------------------------------- */ void Domain::init() { // set box_change flags if box size/shape/sub-domains ever change // due to shrink-wrapping or fixes that change box size/shape/sub-domains box_change_size = box_change_shape = box_change_domain = 0; if (nonperiodic == 2) box_change_size = 1; for (int i = 0; i < modify->nfix; i++) { if (modify->fix[i]->box_change_size) box_change_size = 1; if (modify->fix[i]->box_change_shape) box_change_shape = 1; if (modify->fix[i]->box_change_domain) box_change_domain = 1; } box_change = 0; if (box_change_size || box_change_shape || box_change_domain) box_change = 1; // check for fix deform deform_flag = deform_vremap = deform_groupbit = 0; for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { deform_flag = 1; if (((FixDeform *) modify->fix[i])->remapflag == V_REMAP) { deform_vremap = 1; deform_groupbit = modify->fix[i]->groupbit; } } // region inits for (int i = 0; i < nregion; i++) regions[i]->init(); } /* ---------------------------------------------------------------------- set initial global box assumes boxlo/hi and triclinic tilts are already set expandflag = 1 if need to expand box in shrink-wrapped dims not invoked by read_restart since box is already expanded if don't prevent further expansion, restarted triclinic box with unchanged tilt factors can become a box with atoms outside the box ------------------------------------------------------------------------- */ void Domain::set_initial_box(int expandflag) { // error checks for orthogonal and triclinic domains if (boxlo[0] >= boxhi[0] || boxlo[1] >= boxhi[1] || boxlo[2] >= boxhi[2]) error->one(FLERR,"Box bounds are invalid"); if (domain->dimension == 2 && (xz != 0.0 || yz != 0.0)) error->all(FLERR,"Cannot skew triclinic box in z for 2d simulation"); // error check or warning on triclinic tilt factors if (triclinic) { if ((fabs(xy/(boxhi[0]-boxlo[0])) > 0.5 && xperiodic) || (fabs(xz/(boxhi[0]-boxlo[0])) > 0.5 && xperiodic) || (fabs(yz/(boxhi[1]-boxlo[1])) > 0.5 && yperiodic)) { if (tiltsmall) error->all(FLERR,"Triclinic box skew is too large"); else if (comm->me == 0) error->warning(FLERR,"Triclinic box skew is large"); } } // set small based on box size and SMALL // this works for any unit system small[0] = SMALL * (boxhi[0] - boxlo[0]); small[1] = SMALL * (boxhi[1] - boxlo[1]); small[2] = SMALL * (boxhi[2] - boxlo[2]); // if expandflag, adjust box lo/hi for shrink-wrapped dims if (!expandflag) return; if (boundary[0][0] == 2) boxlo[0] -= small[0]; else if (boundary[0][0] == 3) minxlo = boxlo[0]; if (boundary[0][1] == 2) boxhi[0] += small[0]; else if (boundary[0][1] == 3) minxhi = boxhi[0]; if (boundary[1][0] == 2) boxlo[1] -= small[1]; else if (boundary[1][0] == 3) minylo = boxlo[1]; if (boundary[1][1] == 2) boxhi[1] += small[1]; else if (boundary[1][1] == 3) minyhi = boxhi[1]; if (boundary[2][0] == 2) boxlo[2] -= small[2]; else if (boundary[2][0] == 3) minzlo = boxlo[2]; if (boundary[2][1] == 2) boxhi[2] += small[2]; else if (boundary[2][1] == 3) minzhi = boxhi[2]; } /* ---------------------------------------------------------------------- set global box params assumes boxlo/hi and triclinic tilts are already set ------------------------------------------------------------------------- */ void Domain::set_global_box() { prd[0] = xprd = boxhi[0] - boxlo[0]; prd[1] = yprd = boxhi[1] - boxlo[1]; prd[2] = zprd = boxhi[2] - boxlo[2]; h[0] = xprd; h[1] = yprd; h[2] = zprd; h_inv[0] = 1.0/h[0]; h_inv[1] = 1.0/h[1]; h_inv[2] = 1.0/h[2]; prd_half[0] = xprd_half = 0.5*xprd; prd_half[1] = yprd_half = 0.5*yprd; prd_half[2] = zprd_half = 0.5*zprd; if (triclinic) { h[3] = yz; h[4] = xz; h[5] = xy; h_inv[3] = -h[3] / (h[1]*h[2]); h_inv[4] = (h[3]*h[5] - h[1]*h[4]) / (h[0]*h[1]*h[2]); h_inv[5] = -h[5] / (h[0]*h[1]); boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy); boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz); boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz); boxlo_bound[2] = boxlo[2]; boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy); boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz); boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz); boxhi_bound[2] = boxhi[2]; } } /* ---------------------------------------------------------------------- set lamda box params assumes global box is defined and proc assignment has been made uses comm->xyz_split or comm->mysplit to define subbox boundaries in consistent manner ------------------------------------------------------------------------- */ void Domain::set_lamda_box() { if (comm->layout != LAYOUT_TILED) { int *myloc = comm->myloc; double *xsplit = comm->xsplit; double *ysplit = comm->ysplit; double *zsplit = comm->zsplit; sublo_lamda[0] = xsplit[myloc[0]]; subhi_lamda[0] = xsplit[myloc[0]+1]; sublo_lamda[1] = ysplit[myloc[1]]; subhi_lamda[1] = ysplit[myloc[1]+1]; sublo_lamda[2] = zsplit[myloc[2]]; subhi_lamda[2] = zsplit[myloc[2]+1]; } else { double (*mysplit)[2] = comm->mysplit; sublo_lamda[0] = mysplit[0][0]; subhi_lamda[0] = mysplit[0][1]; sublo_lamda[1] = mysplit[1][0]; subhi_lamda[1] = mysplit[1][1]; sublo_lamda[2] = mysplit[2][0]; subhi_lamda[2] = mysplit[2][1]; } } /* ---------------------------------------------------------------------- set local subbox params for orthogonal boxes assumes global box is defined and proc assignment has been made uses comm->xyz_split or comm->mysplit to define subbox boundaries in consistent manner insure subhi[max] = boxhi ------------------------------------------------------------------------- */ void Domain::set_local_box() { if (triclinic) return; if (comm->layout != LAYOUT_TILED) { int *myloc = comm->myloc; int *procgrid = comm->procgrid; double *xsplit = comm->xsplit; double *ysplit = comm->ysplit; double *zsplit = comm->zsplit; sublo[0] = boxlo[0] + xprd*xsplit[myloc[0]]; if (myloc[0] < procgrid[0]-1) subhi[0] = boxlo[0] + xprd*xsplit[myloc[0]+1]; else subhi[0] = boxhi[0]; sublo[1] = boxlo[1] + yprd*ysplit[myloc[1]]; if (myloc[1] < procgrid[1]-1) subhi[1] = boxlo[1] + yprd*ysplit[myloc[1]+1]; else subhi[1] = boxhi[1]; sublo[2] = boxlo[2] + zprd*zsplit[myloc[2]]; if (myloc[2] < procgrid[2]-1) subhi[2] = boxlo[2] + zprd*zsplit[myloc[2]+1]; else subhi[2] = boxhi[2]; } else { double (*mysplit)[2] = comm->mysplit; sublo[0] = boxlo[0] + xprd*mysplit[0][0]; if (mysplit[0][1] < 1.0) subhi[0] = boxlo[0] + xprd*mysplit[0][1]; else subhi[0] = boxhi[0]; sublo[1] = boxlo[1] + yprd*mysplit[1][0]; if (mysplit[1][1] < 1.0) subhi[1] = boxlo[1] + yprd*mysplit[1][1]; else subhi[1] = boxhi[1]; sublo[2] = boxlo[2] + zprd*mysplit[2][0]; if (mysplit[2][1] < 1.0) subhi[2] = boxlo[2] + zprd*mysplit[2][1]; else subhi[2] = boxhi[2]; } } /* ---------------------------------------------------------------------- reset global & local boxes due to global box boundary changes if shrink-wrapped, determine atom extent and reset boxlo/hi for triclinic, atoms must be in lamda coords (0-1) before reset_box is called ------------------------------------------------------------------------- */ void Domain::reset_box() { // perform shrink-wrapping // compute extent of atoms on this proc // for triclinic, this is done in lamda space if (nonperiodic == 2) { double extent[3][2],all[3][2]; extent[2][0] = extent[1][0] = extent[0][0] = BIG; extent[2][1] = extent[1][1] = extent[0][1] = -BIG; double **x = atom->x; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { extent[0][0] = MIN(extent[0][0],x[i][0]); extent[0][1] = MAX(extent[0][1],x[i][0]); extent[1][0] = MIN(extent[1][0],x[i][1]); extent[1][1] = MAX(extent[1][1],x[i][1]); extent[2][0] = MIN(extent[2][0],x[i][2]); extent[2][1] = MAX(extent[2][1],x[i][2]); } // compute extent across all procs // flip sign of MIN to do it in one Allreduce MAX extent[0][0] = -extent[0][0]; extent[1][0] = -extent[1][0]; extent[2][0] = -extent[2][0]; MPI_Allreduce(extent,all,6,MPI_DOUBLE,MPI_MAX,world); // for triclinic, convert back to box coords before changing box if (triclinic) lamda2x(atom->nlocal); // in shrink-wrapped dims, set box by atom extent // if minimum set, enforce min box size settings // for triclinic, convert lamda extent to box coords, then set box lo/hi // decided NOT to do the next comment - don't want to sneakily change tilt // for triclinic, adjust tilt factors if 2nd dim is shrink-wrapped, // so that displacement in 1st dim stays the same if (triclinic == 0) { if (xperiodic == 0) { if (boundary[0][0] == 2) boxlo[0] = -all[0][0] - small[0]; else if (boundary[0][0] == 3) boxlo[0] = MIN(-all[0][0]-small[0],minxlo); if (boundary[0][1] == 2) boxhi[0] = all[0][1] + small[0]; else if (boundary[0][1] == 3) boxhi[0] = MAX(all[0][1]+small[0],minxhi); if (boxlo[0] > boxhi[0]) error->all(FLERR,"Illegal simulation box"); } if (yperiodic == 0) { if (boundary[1][0] == 2) boxlo[1] = -all[1][0] - small[1]; else if (boundary[1][0] == 3) boxlo[1] = MIN(-all[1][0]-small[1],minylo); if (boundary[1][1] == 2) boxhi[1] = all[1][1] + small[1]; else if (boundary[1][1] == 3) boxhi[1] = MAX(all[1][1]+small[1],minyhi); if (boxlo[1] > boxhi[1]) error->all(FLERR,"Illegal simulation box"); } if (zperiodic == 0) { if (boundary[2][0] == 2) boxlo[2] = -all[2][0] - small[2]; else if (boundary[2][0] == 3) boxlo[2] = MIN(-all[2][0]-small[2],minzlo); if (boundary[2][1] == 2) boxhi[2] = all[2][1] + small[2]; else if (boundary[2][1] == 3) boxhi[2] = MAX(all[2][1]+small[2],minzhi); if (boxlo[2] > boxhi[2]) error->all(FLERR,"Illegal simulation box"); } } else { double lo[3],hi[3]; if (xperiodic == 0) { lo[0] = -all[0][0]; lo[1] = 0.0; lo[2] = 0.0; lamda2x(lo,lo); hi[0] = all[0][1]; hi[1] = 0.0; hi[2] = 0.0; lamda2x(hi,hi); if (boundary[0][0] == 2) boxlo[0] = lo[0] - small[0]; else if (boundary[0][0] == 3) boxlo[0] = MIN(lo[0]-small[0],minxlo); if (boundary[0][1] == 2) boxhi[0] = hi[0] + small[0]; else if (boundary[0][1] == 3) boxhi[0] = MAX(hi[0]+small[0],minxhi); if (boxlo[0] > boxhi[0]) error->all(FLERR,"Illegal simulation box"); } if (yperiodic == 0) { lo[0] = 0.0; lo[1] = -all[1][0]; lo[2] = 0.0; lamda2x(lo,lo); hi[0] = 0.0; hi[1] = all[1][1]; hi[2] = 0.0; lamda2x(hi,hi); if (boundary[1][0] == 2) boxlo[1] = lo[1] - small[1]; else if (boundary[1][0] == 3) boxlo[1] = MIN(lo[1]-small[1],minylo); if (boundary[1][1] == 2) boxhi[1] = hi[1] + small[1]; else if (boundary[1][1] == 3) boxhi[1] = MAX(hi[1]+small[1],minyhi); if (boxlo[1] > boxhi[1]) error->all(FLERR,"Illegal simulation box"); //xy *= (boxhi[1]-boxlo[1]) / yprd; } if (zperiodic == 0) { lo[0] = 0.0; lo[1] = 0.0; lo[2] = -all[2][0]; lamda2x(lo,lo); hi[0] = 0.0; hi[1] = 0.0; hi[2] = all[2][1]; lamda2x(hi,hi); if (boundary[2][0] == 2) boxlo[2] = lo[2] - small[2]; else if (boundary[2][0] == 3) boxlo[2] = MIN(lo[2]-small[2],minzlo); if (boundary[2][1] == 2) boxhi[2] = hi[2] + small[2]; else if (boundary[2][1] == 3) boxhi[2] = MAX(hi[2]+small[2],minzhi); if (boxlo[2] > boxhi[2]) error->all(FLERR,"Illegal simulation box"); //xz *= (boxhi[2]-boxlo[2]) / xprd; //yz *= (boxhi[2]-boxlo[2]) / yprd; } } } // reset box whether shrink-wrapping or not set_global_box(); set_local_box(); // if shrink-wrapped & kspace is defined (i.e. using MSM), call setup() // also call init() (to test for compatibility) ? if (nonperiodic == 2 && force->kspace) { //force->kspace->init(); force->kspace->setup(); } // if shrink-wrapped & triclinic, re-convert to lamda coords for new box // re-invoke pbc() b/c x2lamda result can be outside [0,1] due to roundoff if (nonperiodic == 2 && triclinic) { x2lamda(atom->nlocal); pbc(); } } /* ---------------------------------------------------------------------- enforce PBC and modify box image flags for each atom called every reneighboring and by other commands that change atoms resulting coord must satisfy lo <= coord < hi MAX is important since coord - prd < lo can happen when coord = hi if fix deform, remap velocity of fix group atoms by box edge velocities for triclinic, atoms must be in lamda coords (0-1) before pbc is called image = 10 or 20 bits for each dimension depending on sizeof(imageint) increment/decrement in wrap-around fashion ------------------------------------------------------------------------- */ void Domain::pbc() { int i; imageint idim,otherdims; double *lo,*hi,*period; int nlocal = atom->nlocal; double **x = atom->x; double **v = atom->v; int *mask = atom->mask; imageint *image = atom->image; if (triclinic == 0) { lo = boxlo; hi = boxhi; period = prd; } else { lo = boxlo_lamda; hi = boxhi_lamda; period = prd_lamda; } for (i = 0; i < nlocal; i++) { if (xperiodic) { if (x[i][0] < lo[0]) { x[i][0] += period[0]; if (deform_vremap && mask[i] & deform_groupbit) v[i][0] += h_rate[0]; idim = image[i] & IMGMASK; otherdims = image[i] ^ idim; idim--; idim &= IMGMASK; image[i] = otherdims | idim; } if (x[i][0] >= hi[0]) { x[i][0] -= period[0]; x[i][0] = MAX(x[i][0],lo[0]); if (deform_vremap && mask[i] & deform_groupbit) v[i][0] -= h_rate[0]; idim = image[i] & IMGMASK; otherdims = image[i] ^ idim; idim++; idim &= IMGMASK; image[i] = otherdims | idim; } } if (yperiodic) { if (x[i][1] < lo[1]) { x[i][1] += period[1]; if (deform_vremap && mask[i] & deform_groupbit) { v[i][0] += h_rate[5]; v[i][1] += h_rate[1]; } idim = (image[i] >> IMGBITS) & IMGMASK; otherdims = image[i] ^ (idim << IMGBITS); idim--; idim &= IMGMASK; image[i] = otherdims | (idim << IMGBITS); } if (x[i][1] >= hi[1]) { x[i][1] -= period[1]; x[i][1] = MAX(x[i][1],lo[1]); if (deform_vremap && mask[i] & deform_groupbit) { v[i][0] -= h_rate[5]; v[i][1] -= h_rate[1]; } idim = (image[i] >> IMGBITS) & IMGMASK; otherdims = image[i] ^ (idim << IMGBITS); idim++; idim &= IMGMASK; image[i] = otherdims | (idim << IMGBITS); } } if (zperiodic) { if (x[i][2] < lo[2]) { x[i][2] += period[2]; if (deform_vremap && mask[i] & deform_groupbit) { v[i][0] += h_rate[4]; v[i][1] += h_rate[3]; v[i][2] += h_rate[2]; } idim = image[i] >> IMG2BITS; otherdims = image[i] ^ (idim << IMG2BITS); idim--; idim &= IMGMASK; image[i] = otherdims | (idim << IMG2BITS); } if (x[i][2] >= hi[2]) { x[i][2] -= period[2]; x[i][2] = MAX(x[i][2],lo[2]); if (deform_vremap && mask[i] & deform_groupbit) { v[i][0] -= h_rate[4]; v[i][1] -= h_rate[3]; v[i][2] -= h_rate[2]; } idim = image[i] >> IMG2BITS; otherdims = image[i] ^ (idim << IMG2BITS); idim++; idim &= IMGMASK; image[i] = otherdims | (idim << IMG2BITS); } } } } /* ---------------------------------------------------------------------- check that point is inside box boundaries, in [lo,hi) sense return 1 if true, 0 if false ------------------------------------------------------------------------- */ int Domain::inside(double* x) { double *lo,*hi; double delta[3]; if (triclinic == 0) { lo = boxlo; hi = boxhi; } else { lo = boxlo_lamda; hi = boxhi_lamda; delta[0] = x[0] - boxlo[0]; delta[1] = x[1] - boxlo[1]; delta[2] = x[2] - boxlo[2]; x[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2]; x[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2]; x[2] = h_inv[2]*delta[2]; } if (x[0] < lo[0] || x[0] >= hi[0] || x[1] < lo[1] || x[1] >= hi[1] || x[2] < lo[2] || x[2] >= hi[2]) return 0; else return 1; } /* ---------------------------------------------------------------------- check that point is inside nonperiodic boundaries, in [lo,hi) sense return 1 if true, 0 if false ------------------------------------------------------------------------- */ int Domain::inside_nonperiodic(double* x) { double *lo,*hi; double delta[3]; if (xperiodic && yperiodic && zperiodic) return 1; if (triclinic == 0) { lo = boxlo; hi = boxhi; } else { lo = boxlo_lamda; hi = boxhi_lamda; delta[0] = x[0] - boxlo[0]; delta[1] = x[1] - boxlo[1]; delta[2] = x[2] - boxlo[2]; x[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2]; x[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2]; x[2] = h_inv[2]*delta[2]; } if (!xperiodic && (x[0] < lo[0] || x[0] >= hi[0])) return 0; if (!yperiodic && (x[1] < lo[1] || x[1] >= hi[1])) return 0; if (!zperiodic && (x[2] < lo[2] || x[2] >= hi[2])) return 0; return 1; } /* ---------------------------------------------------------------------- warn if image flags of any bonded atoms are inconsistent could be a problem when using replicate or fix rigid ------------------------------------------------------------------------- */ void Domain::image_check() { int i,j,k,n,imol,iatom; tagint tagprev; // only need to check if system is molecular and some dimension is periodic // if running verlet/split, don't check on KSpace partition since // it has no ghost atoms and thus bond partners won't exist if (!atom->molecular) return; if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return; if (strcmp(update->integrate_style,"verlet/split") == 0 && universe->iworld != 0) return; // communicate unwrapped position of owned atoms to ghost atoms double **unwrap; memory->create(unwrap,atom->nmax,3,"domain:unwrap"); double **x = atom->x; imageint *image = atom->image; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) unmap(x[i],image[i],unwrap[i]); comm->forward_comm_array(3,unwrap); // compute unwrapped extent of each bond // flag if any bond component is longer than 1/2 of periodic box length // flag if any bond component is longer than non-periodic box length // which means image flags in that dimension were different int molecular = atom->molecular; int *num_bond = atom->num_bond; tagint **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; tagint *tag = atom->tag; int *molindex = atom->molindex; int *molatom = atom->molatom; Molecule **onemols = atom->avec->onemols; double delx,dely,delz; int lostbond = output->thermo->lostbond; int nmissing = 0; int flag = 0; for (i = 0; i < nlocal; i++) { if (molecular == 1) n = num_bond[i]; else { if (molindex[i] < 0) continue; imol = molindex[i]; iatom = molatom[i]; n = onemols[imol]->num_bond[iatom]; } for (j = 0; j < n; j++) { if (molecular == 1) { if (bond_type[i][j] <= 0) continue; k = atom->map(bond_atom[i][j]); } else { if (onemols[imol]->bond_type[iatom][j] < 0) continue; tagprev = tag[i] - iatom - 1; k = atom->map(onemols[imol]->bond_atom[iatom][j]+tagprev); } if (k == -1) { nmissing++; if (lostbond == ERROR) error->one(FLERR,"Bond atom missing in image check"); continue; } delx = unwrap[i][0] - unwrap[k][0]; dely = unwrap[i][1] - unwrap[k][1]; delz = unwrap[i][2] - unwrap[k][2]; if (xperiodic && delx > xprd_half) flag = 1; if (xperiodic && dely > yprd_half) flag = 1; if (dimension == 3 && zperiodic && delz > zprd_half) flag = 1; if (!xperiodic && delx > xprd) flag = 1; if (!yperiodic && dely > yprd) flag = 1; if (dimension == 3 && !zperiodic && delz > zprd) flag = 1; } } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world); if (flagall && comm->me == 0) error->warning(FLERR,"Inconsistent image flags"); if (lostbond == WARN) { int all; MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world); if (all && comm->me == 0) error->warning(FLERR,"Bond atom missing in image check"); } memory->destroy(unwrap); } /* ---------------------------------------------------------------------- warn if end atoms in any bonded interaction are further apart than half a periodic box length could cause problems when bonded neighbor list is built since closest_image() could return wrong image ------------------------------------------------------------------------- */ void Domain::box_too_small_check() { int i,j,k,n,imol,iatom; tagint tagprev; // only need to check if system is molecular and some dimension is periodic // if running verlet/split, don't check on KSpace partition since // it has no ghost atoms and thus bond partners won't exist if (!atom->molecular) return; if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return; if (strcmp(update->integrate_style,"verlet/split") == 0 && universe->iworld != 0) return; // maxbondall = longest current bond length // if periodic box dim is tiny (less than 2 * bond-length), // minimum_image() itself may compute bad bond lengths // in this case, image_check() should warn, // assuming 2 atoms have consistent image flags int molecular = atom->molecular; double **x = atom->x; int *num_bond = atom->num_bond; tagint **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; tagint *tag = atom->tag; int *molindex = atom->molindex; int *molatom = atom->molatom; Molecule **onemols = atom->avec->onemols; int nlocal = atom->nlocal; double delx,dely,delz,rsq; double maxbondme = 0.0; int lostbond = output->thermo->lostbond; int nmissing = 0; for (i = 0; i < nlocal; i++) { if (molecular == 1) n = num_bond[i]; else { if (molindex[i] < 0) continue; imol = molindex[i]; iatom = molatom[i]; n = onemols[imol]->num_bond[iatom]; } for (j = 0; j < n; j++) { if (molecular == 1) { if (bond_type[i][j] <= 0) continue; k = atom->map(bond_atom[i][j]); } else { if (onemols[imol]->bond_type[iatom][j] < 0) continue; tagprev = tag[i] - iatom - 1; k = atom->map(onemols[imol]->bond_atom[iatom][j]+tagprev); } if (k == -1) { nmissing++; if (lostbond == ERROR) error->one(FLERR,"Bond atom missing in box size check"); continue; } delx = x[i][0] - x[k][0]; dely = x[i][1] - x[k][1]; delz = x[i][2] - x[k][2]; minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; maxbondme = MAX(maxbondme,rsq); } } if (lostbond == WARN) { int all; MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world); if (all && comm->me == 0) error->warning(FLERR,"Bond atom missing in box size check"); } double maxbondall; MPI_Allreduce(&maxbondme,&maxbondall,1,MPI_DOUBLE,MPI_MAX,world); maxbondall = sqrt(maxbondall); // maxdelta = furthest apart 2 atoms in a bonded interaction can be // include BONDSTRETCH factor to account for dynamics double maxdelta = maxbondall * BONDSTRETCH; if (atom->nangles) maxdelta = 2.0 * maxbondall * BONDSTRETCH; if (atom->ndihedrals) maxdelta = 3.0 * maxbondall * BONDSTRETCH; // warn if maxdelta > than half any periodic box length // since atoms in the interaction could rotate into that dimension int flag = 0; if (xperiodic && maxdelta > xprd_half) flag = 1; if (yperiodic && maxdelta > yprd_half) flag = 1; if (dimension == 3 && zperiodic && maxdelta > zprd_half) flag = 1; int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world); if (flagall && comm->me == 0) error->warning(FLERR, "Bond/angle/dihedral extent > half of periodic box length"); } /* ---------------------------------------------------------------------- check warn if any proc's subbox is smaller than thresh since may lead to lost atoms in comm->exchange() current callers set thresh = neighbor skin ------------------------------------------------------------------------- */ void Domain::subbox_too_small_check(double thresh) { int flag = 0; if (!triclinic) { if (subhi[0]-sublo[0] < thresh || subhi[1]-sublo[1] < thresh) flag = 1; if (dimension == 3 && subhi[2]-sublo[2] < thresh) flag = 1; } else { double delta = subhi_lamda[0] - sublo_lamda[0]; if (delta*prd[0] < thresh) flag = 1; delta = subhi_lamda[1] - sublo_lamda[1]; if (delta*prd[1] < thresh) flag = 1; if (dimension == 3) { delta = subhi_lamda[2] - sublo_lamda[2]; if (delta*prd[2] < thresh) flag = 1; } } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall && comm->me == 0) error->warning(FLERR,"Proc sub-domain size < neighbor skin, " "could lead to lost atoms"); } /* ---------------------------------------------------------------------- minimum image convention use 1/2 of box size as test for triclinic, also add/subtract tilt factors in other dims as needed ------------------------------------------------------------------------- */ void Domain::minimum_image(double &dx, double &dy, double &dz) { if (triclinic == 0) { if (xperiodic) { if (fabs(dx) > xprd_half) { if (dx < 0.0) dx += xprd; else dx -= xprd; } } if (yperiodic) { if (fabs(dy) > yprd_half) { if (dy < 0.0) dy += yprd; else dy -= yprd; } } if (zperiodic) { if (fabs(dz) > zprd_half) { if (dz < 0.0) dz += zprd; else dz -= zprd; } } } else { if (zperiodic) { if (fabs(dz) > zprd_half) { if (dz < 0.0) { dz += zprd; dy += yz; dx += xz; } else { dz -= zprd; dy -= yz; dx -= xz; } } } if (yperiodic) { if (fabs(dy) > yprd_half) { if (dy < 0.0) { dy += yprd; dx += xy; } else { dy -= yprd; dx -= xy; } } } if (xperiodic) { if (fabs(dx) > xprd_half) { if (dx < 0.0) dx += xprd; else dx -= xprd; } } } } /* ---------------------------------------------------------------------- minimum image convention use 1/2 of box size as test for triclinic, also add/subtract tilt factors in other dims as needed ------------------------------------------------------------------------- */ void Domain::minimum_image(double *delta) { if (triclinic == 0) { if (xperiodic) { if (fabs(delta[0]) > xprd_half) { if (delta[0] < 0.0) delta[0] += xprd; else delta[0] -= xprd; } } if (yperiodic) { if (fabs(delta[1]) > yprd_half) { if (delta[1] < 0.0) delta[1] += yprd; else delta[1] -= yprd; } } if (zperiodic) { if (fabs(delta[2]) > zprd_half) { if (delta[2] < 0.0) delta[2] += zprd; else delta[2] -= zprd; } } } else { if (zperiodic) { if (fabs(delta[2]) > zprd_half) { if (delta[2] < 0.0) { delta[2] += zprd; delta[1] += yz; delta[0] += xz; } else { delta[2] -= zprd; delta[1] -= yz; delta[0] -= xz; } } } if (yperiodic) { if (fabs(delta[1]) > yprd_half) { if (delta[1] < 0.0) { delta[1] += yprd; delta[0] += xy; } else { delta[1] -= yprd; delta[0] -= xy; } } } if (xperiodic) { if (fabs(delta[0]) > xprd_half) { if (delta[0] < 0.0) delta[0] += xprd; else delta[0] -= xprd; } } } } /* ---------------------------------------------------------------------- return local index of atom J or any of its images that is closest to atom I if J is not a valid index like -1, just return it ------------------------------------------------------------------------- */ int Domain::closest_image(int i, int j) { if (j < 0) return j; int *sametag = atom->sametag; double **x = atom->x; double *xi = x[i]; int closest = j; double delx = xi[0] - x[j][0]; double dely = xi[1] - x[j][1]; double delz = xi[2] - x[j][2]; double rsqmin = delx*delx + dely*dely + delz*delz; double rsq; while (sametag[j] >= 0) { j = sametag[j]; delx = xi[0] - x[j][0]; dely = xi[1] - x[j][1]; delz = xi[2] - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < rsqmin) { rsqmin = rsq; closest = j; } } return closest; } /* ---------------------------------------------------------------------- find and return Xj image = periodic image of Xj that is closest to Xi for triclinic, add/subtract tilt factors in other dims as needed ------------------------------------------------------------------------- */ void Domain::closest_image(const double * const xi, const double * const xj, double * const xjimage) { double dx = xj[0] - xi[0]; double dy = xj[1] - xi[1]; double dz = xj[2] - xi[2]; if (triclinic == 0) { if (xperiodic) { if (dx < 0.0) { while (dx < 0.0) dx += xprd; if (dx > xprd_half) dx -= xprd; } else { while (dx > 0.0) dx -= xprd; if (dx < -xprd_half) dx += xprd; } } if (yperiodic) { if (dy < 0.0) { while (dy < 0.0) dy += yprd; if (dy > yprd_half) dy -= yprd; } else { while (dy > 0.0) dy -= yprd; if (dy < -yprd_half) dy += yprd; } } if (zperiodic) { if (dz < 0.0) { while (dz < 0.0) dz += zprd; if (dz > zprd_half) dz -= zprd; } else { while (dz > 0.0) dz -= zprd; if (dz < -zprd_half) dz += zprd; } } } else { if (zperiodic) { if (dz < 0.0) { while (dz < 0.0) { dz += zprd; dy += yz; dx += xz; } if (dz > zprd_half) { dz -= zprd; dy -= yz; dx -= xz; } } else { while (dz > 0.0) { dz -= zprd; dy -= yz; dx -= xz; } if (dz < -zprd_half) { dz += zprd; dy += yz; dx += xz; } } } if (yperiodic) { if (dy < 0.0) { while (dy < 0.0) { dy += yprd; dx += xy; } if (dy > yprd_half) { dy -= yprd; dx -= xy; } } else { while (dy > 0.0) { dy -= yprd; dx -= xy; } if (dy < -yprd_half) { dy += yprd; dx += xy; } } } if (xperiodic) { if (dx < 0.0) { while (dx < 0.0) dx += xprd; if (dx > xprd_half) dx -= xprd; } else { while (dx > 0.0) dx -= xprd; if (dx < -xprd_half) dx += xprd; } } } xjimage[0] = xi[0] + dx; xjimage[1] = xi[1] + dy; xjimage[2] = xi[2] + dz; } /* ---------------------------------------------------------------------- remap the point into the periodic box no matter how far away adjust 3 image flags encoded in image accordingly resulting coord must satisfy lo <= coord < hi MAX is important since coord - prd < lo can happen when coord = hi for triclinic, point is converted to lamda coords (0-1) before doing remap image = 10 bits for each dimension increment/decrement in wrap-around fashion ------------------------------------------------------------------------- */ void Domain::remap(double *x, imageint &image) { double *lo,*hi,*period,*coord; double lamda[3]; imageint idim,otherdims; if (triclinic == 0) { lo = boxlo; hi = boxhi; period = prd; coord = x; } else { lo = boxlo_lamda; hi = boxhi_lamda; period = prd_lamda; x2lamda(x,lamda); coord = lamda; } if (xperiodic) { while (coord[0] < lo[0]) { coord[0] += period[0]; idim = image & IMGMASK; otherdims = image ^ idim; idim--; idim &= IMGMASK; image = otherdims | idim; } while (coord[0] >= hi[0]) { coord[0] -= period[0]; idim = image & IMGMASK; otherdims = image ^ idim; idim++; idim &= IMGMASK; image = otherdims | idim; } coord[0] = MAX(coord[0],lo[0]); } if (yperiodic) { while (coord[1] < lo[1]) { coord[1] += period[1]; idim = (image >> IMGBITS) & IMGMASK; otherdims = image ^ (idim << IMGBITS); idim--; idim &= IMGMASK; image = otherdims | (idim << IMGBITS); } while (coord[1] >= hi[1]) { coord[1] -= period[1]; idim = (image >> IMGBITS) & IMGMASK; otherdims = image ^ (idim << IMGBITS); idim++; idim &= IMGMASK; image = otherdims | (idim << IMGBITS); } coord[1] = MAX(coord[1],lo[1]); } if (zperiodic) { while (coord[2] < lo[2]) { coord[2] += period[2]; idim = image >> IMG2BITS; otherdims = image ^ (idim << IMG2BITS); idim--; idim &= IMGMASK; image = otherdims | (idim << IMG2BITS); } while (coord[2] >= hi[2]) { coord[2] -= period[2]; idim = image >> IMG2BITS; otherdims = image ^ (idim << IMG2BITS); idim++; idim &= IMGMASK; image = otherdims | (idim << IMG2BITS); } coord[2] = MAX(coord[2],lo[2]); } if (triclinic) lamda2x(coord,x); } /* ---------------------------------------------------------------------- remap the point into the periodic box no matter how far away no image flag calculation resulting coord must satisfy lo <= coord < hi MAX is important since coord - prd < lo can happen when coord = hi for triclinic, point is converted to lamda coords (0-1) before remap ------------------------------------------------------------------------- */ void Domain::remap(double *x) { double *lo,*hi,*period,*coord; double lamda[3]; if (triclinic == 0) { lo = boxlo; hi = boxhi; period = prd; coord = x; } else { lo = boxlo_lamda; hi = boxhi_lamda; period = prd_lamda; x2lamda(x,lamda); coord = lamda; } if (xperiodic) { while (coord[0] < lo[0]) coord[0] += period[0]; while (coord[0] >= hi[0]) coord[0] -= period[0]; coord[0] = MAX(coord[0],lo[0]); } if (yperiodic) { while (coord[1] < lo[1]) coord[1] += period[1]; while (coord[1] >= hi[1]) coord[1] -= period[1]; coord[1] = MAX(coord[1],lo[1]); } if (zperiodic) { while (coord[2] < lo[2]) coord[2] += period[2]; while (coord[2] >= hi[2]) coord[2] -= period[2]; coord[2] = MAX(coord[2],lo[2]); } if (triclinic) lamda2x(coord,x); } /* ---------------------------------------------------------------------- remap xnew to be within half box length of xold do it directly, not iteratively, in case is far away for triclinic, both points are converted to lamda coords (0-1) before remap ------------------------------------------------------------------------- */ void Domain::remap_near(double *xnew, double *xold) { int n; double *coordnew,*coordold,*period,*half; double lamdanew[3],lamdaold[3]; if (triclinic == 0) { period = prd; half = prd_half; coordnew = xnew; coordold = xold; } else { period = prd_lamda; half = prd_half_lamda; x2lamda(xnew,lamdanew); coordnew = lamdanew; x2lamda(xold,lamdaold); coordold = lamdaold; } // iterative form // if (xperiodic) { // while (coordnew[0]-coordold[0] > half[0]) coordnew[0] -= period[0]; // while (coordold[0]-coordnew[0] > half[0]) coordnew[0] += period[0]; // } if (xperiodic) { if (coordnew[0]-coordold[0] > period[0]) { n = static_cast ((coordnew[0]-coordold[0])/period[0]); coordnew[0] -= n*period[0]; } while (coordnew[0]-coordold[0] > half[0]) coordnew[0] -= period[0]; if (coordold[0]-coordnew[0] > period[0]) { n = static_cast ((coordold[0]-coordnew[0])/period[0]); coordnew[0] += n*period[0]; } while (coordold[0]-coordnew[0] > half[0]) coordnew[0] += period[0]; } if (yperiodic) { if (coordnew[1]-coordold[1] > period[1]) { n = static_cast ((coordnew[1]-coordold[1])/period[1]); coordnew[1] -= n*period[1]; } while (coordnew[1]-coordold[1] > half[1]) coordnew[1] -= period[1]; if (coordold[1]-coordnew[1] > period[1]) { n = static_cast ((coordold[1]-coordnew[1])/period[1]); coordnew[1] += n*period[1]; } while (coordold[1]-coordnew[1] > half[1]) coordnew[1] += period[1]; } if (zperiodic) { if (coordnew[2]-coordold[2] > period[2]) { n = static_cast ((coordnew[2]-coordold[2])/period[2]); coordnew[2] -= n*period[2]; } while (coordnew[2]-coordold[2] > half[2]) coordnew[2] -= period[2]; if (coordold[2]-coordnew[2] > period[2]) { n = static_cast ((coordold[2]-coordnew[2])/period[2]); coordnew[2] += n*period[2]; } while (coordold[2]-coordnew[2] > half[2]) coordnew[2] += period[2]; } if (triclinic) lamda2x(coordnew,xnew); } /* ---------------------------------------------------------------------- unmap the point via image flags x overwritten with result, don't reset image flag for triclinic, use h[] to add in tilt factors in other dims as needed ------------------------------------------------------------------------- */ void Domain::unmap(double *x, imageint image) { int xbox = (image & IMGMASK) - IMGMAX; int ybox = (image >> IMGBITS & IMGMASK) - IMGMAX; int zbox = (image >> IMG2BITS) - IMGMAX; if (triclinic == 0) { x[0] += xbox*xprd; x[1] += ybox*yprd; x[2] += zbox*zprd; } else { x[0] += h[0]*xbox + h[5]*ybox + h[4]*zbox; x[1] += h[1]*ybox + h[3]*zbox; x[2] += h[2]*zbox; } } /* ---------------------------------------------------------------------- unmap the point via image flags result returned in y, don't reset image flag for triclinic, use h[] to add in tilt factors in other dims as needed ------------------------------------------------------------------------- */ void Domain::unmap(double *x, imageint image, double *y) { int xbox = (image & IMGMASK) - IMGMAX; int ybox = (image >> IMGBITS & IMGMASK) - IMGMAX; int zbox = (image >> IMG2BITS) - IMGMAX; if (triclinic == 0) { y[0] = x[0] + xbox*xprd; y[1] = x[1] + ybox*yprd; y[2] = x[2] + zbox*zprd; } else { y[0] = x[0] + h[0]*xbox + h[5]*ybox + h[4]*zbox; y[1] = x[1] + h[1]*ybox + h[3]*zbox; y[2] = x[2] + h[2]*zbox; } } /* ---------------------------------------------------------------------- adjust image flags due to triclinic box flip flip operation is changing box vectors A,B,C to new A',B',C' A' = A (A does not change) B' = B + mA (B shifted by A) C' = C + pB + nA (C shifted by B and/or A) this requires the image flags change from (a,b,c) to (a',b',c') so that x_unwrap for each atom is same before/after x_unwrap_before = xlocal + aA + bB + cC x_unwrap_after = xlocal + a'A' + b'B' + c'C' this requires: c' = c b' = b - cp a' = a - (b-cp)m - cn = a - b'm - cn in other words, for xy flip, change in x flag depends on current y flag this is b/c the xy flip dramatically changes which tiled image of simulation box an unwrapped point maps to ------------------------------------------------------------------------- */ void Domain::image_flip(int m, int n, int p) { imageint *image = atom->image; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { int xbox = (image[i] & IMGMASK) - IMGMAX; int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; int zbox = (image[i] >> IMG2BITS) - IMGMAX; ybox -= p*zbox; xbox -= m*ybox + n*zbox; image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) | (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS); } } /* ---------------------------------------------------------------------- create a lattice ------------------------------------------------------------------------- */ void Domain::set_lattice(int narg, char **arg) { if (lattice) delete lattice; lattice = new Lattice(lmp,narg,arg); } /* ---------------------------------------------------------------------- create a new region ------------------------------------------------------------------------- */ void Domain::add_region(int narg, char **arg) { if (narg < 2) error->all(FLERR,"Illegal region command"); if (strcmp(arg[1],"delete") == 0) { delete_region(narg,arg); return; } if (find_region(arg[0]) >= 0) error->all(FLERR,"Reuse of region ID"); // extend Region list if necessary if (nregion == maxregion) { maxregion += DELTAREGION; regions = (Region **) memory->srealloc(regions,maxregion*sizeof(Region *),"domain:regions"); } // create the Region if (strcmp(arg[1],"none") == 0) error->all(FLERR,"Unknown region style"); #define REGION_CLASS #define RegionStyle(key,Class) \ else if (strcmp(arg[1],#key) == 0) \ regions[nregion] = new Class(lmp,narg,arg); #include "style_region.h" #undef REGION_CLASS else error->all(FLERR,"Unknown region style"); // initialize any region variables via init() // in case region is used between runs, e.g. to print a variable regions[nregion]->init(); nregion++; } /* ---------------------------------------------------------------------- delete a region ------------------------------------------------------------------------- */ void Domain::delete_region(int narg, char **arg) { if (narg != 2) error->all(FLERR,"Illegal region command"); int iregion = find_region(arg[0]); if (iregion == -1) error->all(FLERR,"Delete region ID does not exist"); delete regions[iregion]; regions[iregion] = regions[nregion-1]; nregion--; } /* ---------------------------------------------------------------------- return region index if name matches existing region ID return -1 if no such region ------------------------------------------------------------------------- */ int Domain::find_region(char *name) { for (int iregion = 0; iregion < nregion; iregion++) if (strcmp(name,regions[iregion]->id) == 0) return iregion; return -1; } /* ---------------------------------------------------------------------- (re)set boundary settings flag = 0, called from the input script flag = 1, called from change box command ------------------------------------------------------------------------- */ void Domain::set_boundary(int narg, char **arg, int flag) { if (narg != 3) error->all(FLERR,"Illegal boundary command"); char c; for (int idim = 0; idim < 3; idim++) for (int iside = 0; iside < 2; iside++) { if (iside == 0) c = arg[idim][0]; else if (iside == 1 && strlen(arg[idim]) == 1) c = arg[idim][0]; else c = arg[idim][1]; if (c == 'p') boundary[idim][iside] = 0; else if (c == 'f') boundary[idim][iside] = 1; else if (c == 's') boundary[idim][iside] = 2; else if (c == 'm') boundary[idim][iside] = 3; else { if (flag == 0) error->all(FLERR,"Illegal boundary command"); if (flag == 1) error->all(FLERR,"Illegal change_box command"); } } for (int idim = 0; idim < 3; idim++) if ((boundary[idim][0] == 0 && boundary[idim][1]) || (boundary[idim][0] && boundary[idim][1] == 0)) error->all(FLERR,"Both sides of boundary must be periodic"); if (boundary[0][0] == 0) xperiodic = 1; else xperiodic = 0; if (boundary[1][0] == 0) yperiodic = 1; else yperiodic = 0; if (boundary[2][0] == 0) zperiodic = 1; else zperiodic = 0; periodicity[0] = xperiodic; periodicity[1] = yperiodic; periodicity[2] = zperiodic; nonperiodic = 0; if (xperiodic == 0 || yperiodic == 0 || zperiodic == 0) { nonperiodic = 1; if (boundary[0][0] >= 2 || boundary[0][1] >= 2 || boundary[1][0] >= 2 || boundary[1][1] >= 2 || boundary[2][0] >= 2 || boundary[2][1] >= 2) nonperiodic = 2; } } /* ---------------------------------------------------------------------- set domain attributes ------------------------------------------------------------------------- */ void Domain::set_box(int narg, char **arg) { if (narg < 1) error->all(FLERR,"Illegal box command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"tilt") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal box command"); if (strcmp(arg[iarg+1],"small") == 0) tiltsmall = 1; else if (strcmp(arg[iarg+1],"large") == 0) tiltsmall = 0; else error->all(FLERR,"Illegal box command"); iarg += 2; } else error->all(FLERR,"Illegal box command"); } } /* ---------------------------------------------------------------------- print box info, orthogonal or triclinic ------------------------------------------------------------------------- */ void Domain::print_box(const char *str) { if (comm->me == 0) { if (screen) { if (triclinic == 0) fprintf(screen,"%sorthogonal box = (%g %g %g) to (%g %g %g)\n", str,boxlo[0],boxlo[1],boxlo[2],boxhi[0],boxhi[1],boxhi[2]); else { char *format = (char *) "%striclinic box = (%g %g %g) to (%g %g %g) with tilt (%g %g %g)\n"; fprintf(screen,format, str,boxlo[0],boxlo[1],boxlo[2],boxhi[0],boxhi[1],boxhi[2], xy,xz,yz); } } if (logfile) { if (triclinic == 0) fprintf(logfile,"%sorthogonal box = (%g %g %g) to (%g %g %g)\n", str,boxlo[0],boxlo[1],boxlo[2],boxhi[0],boxhi[1],boxhi[2]); else { char *format = (char *) "%striclinic box = (%g %g %g) to (%g %g %g) with tilt (%g %g %g)\n"; fprintf(logfile,format, str,boxlo[0],boxlo[1],boxlo[2],boxhi[0],boxhi[1],boxhi[2], xy,xz,yz); } } } } /* ---------------------------------------------------------------------- format boundary string for output assume str is 9 chars or more in length ------------------------------------------------------------------------- */ void Domain::boundary_string(char *str) { int m = 0; for (int idim = 0; idim < 3; idim++) { for (int iside = 0; iside < 2; iside++) { if (boundary[idim][iside] == 0) str[m++] = 'p'; else if (boundary[idim][iside] == 1) str[m++] = 'f'; else if (boundary[idim][iside] == 2) str[m++] = 's'; else if (boundary[idim][iside] == 3) str[m++] = 'm'; } str[m++] = ' '; } str[8] = '\0'; } /* ---------------------------------------------------------------------- convert triclinic 0-1 lamda coords to box coords for all N atoms x = H lamda + x0; ------------------------------------------------------------------------- */ void Domain::lamda2x(int n) { double **x = atom->x; for (int i = 0; i < n; i++) { x[i][0] = h[0]*x[i][0] + h[5]*x[i][1] + h[4]*x[i][2] + boxlo[0]; x[i][1] = h[1]*x[i][1] + h[3]*x[i][2] + boxlo[1]; x[i][2] = h[2]*x[i][2] + boxlo[2]; } } /* ---------------------------------------------------------------------- convert box coords to triclinic 0-1 lamda coords for all N atoms lamda = H^-1 (x - x0) ------------------------------------------------------------------------- */ void Domain::x2lamda(int n) { double delta[3]; double **x = atom->x; for (int i = 0; i < n; i++) { delta[0] = x[i][0] - boxlo[0]; delta[1] = x[i][1] - boxlo[1]; delta[2] = x[i][2] - boxlo[2]; x[i][0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2]; x[i][1] = h_inv[1]*delta[1] + h_inv[3]*delta[2]; x[i][2] = h_inv[2]*delta[2]; } } /* ---------------------------------------------------------------------- convert triclinic 0-1 lamda coords to box coords for one atom x = H lamda + x0; lamda and x can point to same 3-vector ------------------------------------------------------------------------- */ void Domain::lamda2x(double *lamda, double *x) { x[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2] + boxlo[0]; x[1] = h[1]*lamda[1] + h[3]*lamda[2] + boxlo[1]; x[2] = h[2]*lamda[2] + boxlo[2]; } /* ---------------------------------------------------------------------- convert box coords to triclinic 0-1 lamda coords for one atom lamda = H^-1 (x - x0) x and lamda can point to same 3-vector ------------------------------------------------------------------------- */ void Domain::x2lamda(double *x, double *lamda) { double delta[3]; delta[0] = x[0] - boxlo[0]; delta[1] = x[1] - boxlo[1]; delta[2] = x[2] - boxlo[2]; lamda[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2]; lamda[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2]; lamda[2] = h_inv[2]*delta[2]; } /* ---------------------------------------------------------------------- convert box coords to triclinic 0-1 lamda coords for one atom use my_boxlo & my_h_inv stored by caller for previous state of box lamda = H^-1 (x - x0) x and lamda can point to same 3-vector ------------------------------------------------------------------------- */ void Domain::x2lamda(double *x, double *lamda, double *my_boxlo, double *my_h_inv) { double delta[3]; delta[0] = x[0] - my_boxlo[0]; delta[1] = x[1] - my_boxlo[1]; delta[2] = x[2] - my_boxlo[2]; lamda[0] = my_h_inv[0]*delta[0] + my_h_inv[5]*delta[1] + my_h_inv[4]*delta[2]; lamda[1] = my_h_inv[1]*delta[1] + my_h_inv[3]*delta[2]; lamda[2] = my_h_inv[2]*delta[2]; } /* ---------------------------------------------------------------------- convert 8 lamda corner pts of lo/hi box to box coords return bboxlo/hi = bounding box around 8 corner pts in box coords ------------------------------------------------------------------------- */ void Domain::bbox(double *lo, double *hi, double *bboxlo, double *bboxhi) { double x[3]; bboxlo[0] = bboxlo[1] = bboxlo[2] = BIG; bboxhi[0] = bboxhi[1] = bboxhi[2] = -BIG; x[0] = lo[0]; x[1] = lo[1]; x[2] = lo[2]; lamda2x(x,x); bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]); bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); x[0] = hi[0]; x[1] = lo[1]; x[2] = lo[2]; lamda2x(x,x); bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]); bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); x[0] = lo[0]; x[1] = hi[1]; x[2] = lo[2]; lamda2x(x,x); bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]); bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); x[0] = hi[0]; x[1] = hi[1]; x[2] = lo[2]; lamda2x(x,x); bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]); bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); x[0] = lo[0]; x[1] = lo[1]; x[2] = hi[2]; lamda2x(x,x); bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]); bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); x[0] = hi[0]; x[1] = lo[1]; x[2] = hi[2]; lamda2x(x,x); bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]); bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); x[0] = lo[0]; x[1] = hi[1]; x[2] = hi[2]; lamda2x(x,x); bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]); bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); x[0] = hi[0]; x[1] = hi[1]; x[2] = hi[2]; lamda2x(x,x); bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]); bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); } /* ---------------------------------------------------------------------- compute 8 corner pts of my triclinic sub-box output is in corners, see ordering in lamda_box_corners ------------------------------------------------------------------------- */ void Domain::box_corners() { lamda_box_corners(boxlo_lamda,boxhi_lamda); } /* ---------------------------------------------------------------------- compute 8 corner pts of my triclinic sub-box output is in corners, see ordering in lamda_box_corners ------------------------------------------------------------------------- */ void Domain::subbox_corners() { lamda_box_corners(sublo_lamda,subhi_lamda); } /* ---------------------------------------------------------------------- compute 8 corner pts of any triclinic box with lo/hi in lamda coords 8 output conners are ordered with x changing fastest, then y, finally z could be more efficient if just coded with xy,yz,xz explicitly ------------------------------------------------------------------------- */ void Domain::lamda_box_corners(double *lo, double *hi) { corners[0][0] = lo[0]; corners[0][1] = lo[1]; corners[0][2] = lo[2]; lamda2x(corners[0],corners[0]); corners[1][0] = hi[0]; corners[1][1] = lo[1]; corners[1][2] = lo[2]; lamda2x(corners[1],corners[1]); corners[2][0] = lo[0]; corners[2][1] = hi[1]; corners[2][2] = lo[2]; lamda2x(corners[2],corners[2]); corners[3][0] = hi[0]; corners[3][1] = hi[1]; corners[3][2] = lo[2]; lamda2x(corners[3],corners[3]); corners[4][0] = lo[0]; corners[4][1] = lo[1]; corners[4][2] = hi[2]; lamda2x(corners[4],corners[4]); corners[5][0] = hi[0]; corners[5][1] = lo[1]; corners[5][2] = hi[2]; lamda2x(corners[5],corners[5]); corners[6][0] = lo[0]; corners[6][1] = hi[1]; corners[6][2] = hi[2]; lamda2x(corners[6],corners[6]); corners[7][0] = hi[0]; corners[7][1] = hi[1]; corners[7][2] = subhi_lamda[2]; lamda2x(corners[7],corners[7]); } diff --git a/src/domain.h b/src/domain.h index 1a3ba0834..65681a513 100644 --- a/src/domain.h +++ b/src/domain.h @@ -1,267 +1,269 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifndef LMP_DOMAIN_H #define LMP_DOMAIN_H #include "math.h" #include "pointers.h" namespace LAMMPS_NS { class Domain : protected Pointers { public: int box_exist; // 0 = not yet created, 1 = exists int dimension; // 2 = 2d, 3 = 3d int nonperiodic; // 0 = periodic in all 3 dims // 1 = periodic or fixed in all 6 // 2 = shrink-wrap in any of 6 int xperiodic,yperiodic,zperiodic; // 0 = non-periodic, 1 = periodic int periodicity[3]; // xyz periodicity as array int boundary[3][2]; // settings for 6 boundaries // 0 = periodic // 1 = fixed non-periodic // 2 = shrink-wrap non-periodic // 3 = shrink-wrap non-per w/ min int triclinic; // 0 = orthog box, 1 = triclinic int tiltsmall; // 1 if limit tilt, else 0 // orthogonal box double xprd,yprd,zprd; // global box dimensions double xprd_half,yprd_half,zprd_half; // half dimensions double prd[3]; // array form of dimensions double prd_half[3]; // array form of half dimensions // triclinic box // xprd,xprd_half,prd,prd_half = // same as if untilted double prd_lamda[3]; // lamda box = (1,1,1) double prd_half_lamda[3]; // lamda half box = (0.5,0.5,0.5) double boxlo[3],boxhi[3]; // orthogonal box global bounds // triclinic box // boxlo/hi = same as if untilted double boxlo_lamda[3],boxhi_lamda[3]; // lamda box = (0,1) double boxlo_bound[3],boxhi_bound[3]; // bounding box of tilted domain double corners[8][3]; // 8 corner points // orthogonal box & triclinic box double minxlo,minxhi; // minimum size of global box double minylo,minyhi; // when shrink-wrapping double minzlo,minzhi; // tri only possible for non-skew dims // orthogonal box double sublo[3],subhi[3]; // sub-box bounds on this proc // triclinic box // sublo/hi = undefined double sublo_lamda[3],subhi_lamda[3]; // bounds of subbox in lamda // triclinic box double xy,xz,yz; // 3 tilt factors double h[6],h_inv[6]; // shape matrix in Voigt notation double h_rate[6],h_ratelo[3]; // rate of box size/shape change int box_change; // 1 if any of next 3 flags are set, else 0 int box_change_size; // 1 if box size changes, 0 if not int box_change_shape; // 1 if box shape changes, 0 if not int box_change_domain; // 1 if proc sub-domains change, 0 if not int deform_flag; // 1 if fix deform exist, else 0 int deform_vremap; // 1 if fix deform remaps v, else 0 int deform_groupbit; // atom group to perform v remap for class Lattice *lattice; // user-defined lattice int nregion; // # of defined Regions int maxregion; // max # list can hold class Region **regions; // list of defined Regions + int copymode; + Domain(class LAMMPS *); virtual ~Domain(); virtual void init(); void set_initial_box(int expandflag=1); virtual void set_global_box(); virtual void set_lamda_box(); virtual void set_local_box(); virtual void reset_box(); virtual void pbc(); void image_check(); void box_too_small_check(); void subbox_too_small_check(double); void minimum_image(double &, double &, double &); void minimum_image(double *); int closest_image(int, int); void closest_image(const double * const, const double * const, double * const); void remap(double *, imageint &); void remap(double *); void remap_near(double *, double *); void unmap(double *, imageint); void unmap(double *, imageint, double *); void image_flip(int, int, int); void set_lattice(int, char **); void add_region(int, char **); void delete_region(int, char **); int find_region(char *); void set_boundary(int, char **, int); void set_box(int, char **); void print_box(const char *); void boundary_string(char *); virtual void lamda2x(int); virtual void x2lamda(int); virtual void lamda2x(double *, double *); virtual void x2lamda(double *, double *); int inside(double *); int inside_nonperiodic(double *); void x2lamda(double *, double *, double *, double *); void bbox(double *, double *, double *, double *); void box_corners(); void subbox_corners(); void lamda_box_corners(double *, double *); // minimum image convention check // return 1 if any distance > 1/2 of box size // indicates a special neighbor is actually not in a bond, // but is a far-away image that should be treated as an unbonded neighbor // inline since called from neighbor build inner loop // inline int minimum_image_check(double dx, double dy, double dz) { if (xperiodic && fabs(dx) > xprd_half) return 1; if (yperiodic && fabs(dy) > yprd_half) return 1; if (zperiodic && fabs(dz) > zprd_half) return 1; return 0; } protected: double small[3]; // fractions of box lengths }; } #endif /* ERROR/WARNING messages: E: Box bounds are invalid The box boundaries specified in the read_data file are invalid. The lo value must be less than the hi value for all 3 dimensions. E: Cannot skew triclinic box in z for 2d simulation Self-explanatory. E: Triclinic box skew is too large The displacement in a skewed direction must be less than half the box length in that dimension. E.g. the xy tilt must be between -half and +half of the x box length. This constraint can be relaxed by using the box tilt command. W: Triclinic box skew is large The displacement in a skewed direction is normally required to be less than half the box length in that dimension. E.g. the xy tilt must be between -half and +half of the x box length. You have relaxed the constraint using the box tilt command, but the warning means that a LAMMPS simulation may be inefficient as a result. E: Illegal simulation box The lower bound of the simulation box is greater than the upper bound. E: Bond atom missing in image check The 2nd atom in a particular bond is missing on this processor. Typically this is because the pairwise cutoff is set too short or the bond has blown apart and an atom is too far away. W: Inconsistent image flags The image flags for a pair on bonded atoms appear to be inconsistent. Inconsistent means that when the coordinates of the two atoms are unwrapped using the image flags, the two atoms are far apart. Specifically they are further apart than half a periodic box length. Or they are more than a box length apart in a non-periodic dimension. This is usually due to the initial data file not having correct image flags for the 2 atoms in a bond that straddles a periodic boundary. They should be different by 1 in that case. This is a warning because inconsistent image flags will not cause problems for dynamics or most LAMMPS simulations. However they can cause problems when such atoms are used with the fix rigid or replicate commands. W: Bond atom missing in image check The 2nd atom in a particular bond is missing on this processor. Typically this is because the pairwise cutoff is set too short or the bond has blown apart and an atom is too far away. E: Bond atom missing in box size check The 2nd atoms needed to compute a particular bond is missing on this processor. Typically this is because the pairwise cutoff is set too short or the bond has blown apart and an atom is too far away. W: Bond atom missing in box size check The 2nd atoms needed to compute a particular bond is missing on this processor. Typically this is because the pairwise cutoff is set too short or the bond has blown apart and an atom is too far away. W: Bond/angle/dihedral extent > half of periodic box length This is a restriction because LAMMPS can be confused about which image of an atom in the bonded interaction is the correct one to use. "Extent" in this context means the maximum end-to-end length of the bond/angle/dihedral. LAMMPS computes this by taking the maximum bond length, multiplying by the number of bonds in the interaction (e.g. 3 for a dihedral) and adding a small amount of stretch. W: Proc sub-domain size < neighbor skin, could lead to lost atoms The decomposition of the physical domain (likely due to load balancing) has led to a processor's sub-domain being smaller than the neighbor skin in one or more dimensions. Since reneighboring is triggered by atoms moving the skin distance, this may lead to lost atoms, if an atom moves all the way across a neighboring processor's sub-domain before reneighboring is triggered. E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Reuse of region ID A region ID cannot be used twice. E: Unknown region style The choice of region style is unknown. E: Delete region ID does not exist Self-explanatory. E: Both sides of boundary must be periodic Cannot specify a boundary as periodic only on the lo or hi side. Must be periodic on both sides. */ diff --git a/src/fix_deform.h b/src/fix_deform.h index 469fa8776..efbecfd63 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -1,136 +1,136 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef FIX_CLASS FixStyle(deform,FixDeform) #else #ifndef LMP_FIX_DEFORM_H #define LMP_FIX_DEFORM_H #include "fix.h" namespace LAMMPS_NS { class FixDeform : public Fix { public: int remapflag; // whether x,v are remapped across PBC int dimflag[6]; // which dims are deformed FixDeform(class LAMMPS *, int, char **); - ~FixDeform(); + virtual ~FixDeform(); int setmask(); void init(); - void pre_exchange(); - void end_of_step(); + virtual void pre_exchange(); + virtual void end_of_step(); double memory_usage(); - private: + protected: int triclinic,scaleflag,flipflag; int flip,flipxy,flipxz,flipyz; double *h_rate,*h_ratelo; int varflag; // 1 if VARIABLE option is used, 0 if not int kspace_flag; // 1 if KSpace invoked, 0 if not int nrigid; // number of rigid fixes int *rfix; // indices of rigid fixes class Irregular *irregular; // for migrating atoms after box flips double TWOPI; struct Set { int style,substyle; double flo,fhi,ftilt; double dlo,dhi,dtilt; double scale,vel,rate; double amplitude,tperiod; double lo_initial,hi_initial; double lo_start,hi_start,lo_stop,hi_stop,lo_target,hi_target; double tilt_initial,tilt_start,tilt_stop,tilt_target,tilt_flip; double tilt_min,tilt_max; double vol_initial,vol_start; int fixed,dynamic1,dynamic2; char *hstr,*hratestr; int hvar,hratevar; }; Set *set; void options(int, char **); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Fix deform tilt factors require triclinic box Cannot deform the tilt factors of a simulation box unless it is a triclinic (non-orthogonal) box. E: Cannot use fix deform on a shrink-wrapped boundary The x, y, z options cannot be applied to shrink-wrapped dimensions. E: Cannot use fix deform tilt on a shrink-wrapped 2nd dim This is because the shrink-wrapping will change the value of the strain implied by the tilt factor. E: Fix deform volume setting is invalid Cannot use volume style unless other dimensions are being controlled. E: More than one fix deform Only one fix deform can be defined at a time. E: Variable name for fix deform does not exist Self-explantory. E: Variable for fix deform is invalid style The variable must be an equal-style variable. E: Final box dimension due to fix deform is < 0.0 Self-explanatory. E: Cannot use fix deform trate on a box with zero tilt The trate style alters the current strain. E: Fix deform cannot use yz variable with xy The yz setting cannot be a variable if xy deformation is also specified. This is because LAMMPS cannot determine if the yz setting will induce a box flip which would be invalid if xy is also changing. E: Fix deform is changing yz too much with xy When both yz and xy are changing, it induces changes in xz if the box must flip from one tilt extreme to another. Thus it is not allowed for yz to grow so much that a flip is induced. */ diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index e9b23a1c9..730b247cd 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -1,2298 +1,2300 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing authors: Mark Stevens (SNL), Aidan Thompson (SNL) ------------------------------------------------------------------------- */ #include "string.h" #include "stdlib.h" #include "math.h" #include "fix_nh.h" #include "math_extra.h" #include "atom.h" #include "force.h" #include "group.h" #include "comm.h" #include "neighbor.h" #include "irregular.h" #include "modify.h" #include "fix_deform.h" #include "compute.h" #include "kspace.h" #include "update.h" #include "respa.h" #include "domain.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; #define DELTAFLIP 0.1 #define TILTMAX 1.5 enum{NOBIAS,BIAS}; enum{NONE,XYZ,XY,YZ,XZ}; enum{ISO,ANISO,TRICLINIC}; /* ---------------------------------------------------------------------- NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion ---------------------------------------------------------------------- */ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command"); restart_global = 1; dynamic_group_allow = 1; time_integrate = 1; scalar_flag = 1; vector_flag = 1; global_freq = 1; extscalar = 1; extvector = 0; // default values pcouple = NONE; drag = 0.0; allremap = 1; id_dilate = NULL; mtchain = mpchain = 3; nc_tchain = nc_pchain = 1; mtk_flag = 1; deviatoric_flag = 0; nreset_h0 = 0; eta_mass_flag = 1; omega_mass_flag = 0; etap_mass_flag = 0; flipflag = 1; // turn on tilt factor scaling, whenever applicable dimension = domain->dimension; scaleyz = scalexz = scalexy = 0; if (domain->yperiodic && domain->xy != 0.0) scalexy = 1; if (domain->zperiodic && dimension == 3) { if (domain->yz != 0.0) scaleyz = 1; if (domain->xz != 0.0) scalexz = 1; } // set fixed-point to default = center of cell fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]); fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]); fixedpoint[2] = 0.5*(domain->boxlo[2]+domain->boxhi[2]); // used by FixNVTSllod to preserve non-default value mtchain_default_flag = 1; tstat_flag = 0; double t_period = 0.0; double p_period[6]; for (int i = 0; i < 6; i++) { p_start[i] = p_stop[i] = p_period[i] = p_target[i] = 0.0; p_flag[i] = 0; } // process keywords int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"temp") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); tstat_flag = 1; t_start = force->numeric(FLERR,arg[iarg+1]); t_target = t_start; t_stop = force->numeric(FLERR,arg[iarg+2]); t_period = force->numeric(FLERR,arg[iarg+3]); if (t_start < 0.0 || t_stop <= 0.0) error->all(FLERR, "Target temperature for fix nvt/npt/nph cannot be 0.0"); iarg += 4; } else if (strcmp(arg[iarg],"iso") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = XYZ; p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = p_flag[1] = p_flag[2] = 1; if (dimension == 2) { p_start[2] = p_stop[2] = p_period[2] = 0.0; p_flag[2] = 0; } iarg += 4; } else if (strcmp(arg[iarg],"aniso") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = NONE; p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = p_flag[1] = p_flag[2] = 1; if (dimension == 2) { p_start[2] = p_stop[2] = p_period[2] = 0.0; p_flag[2] = 0; } iarg += 4; } else if (strcmp(arg[iarg],"tri") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = NONE; scalexy = scalexz = scaleyz = 0; p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = p_flag[1] = p_flag[2] = 1; p_start[3] = p_start[4] = p_start[5] = 0.0; p_stop[3] = p_stop[4] = p_stop[5] = 0.0; p_period[3] = p_period[4] = p_period[5] = force->numeric(FLERR,arg[iarg+3]); p_flag[3] = p_flag[4] = p_flag[5] = 1; if (dimension == 2) { p_start[2] = p_stop[2] = p_period[2] = 0.0; p_flag[2] = 0; p_start[3] = p_stop[3] = p_period[3] = 0.0; p_flag[3] = 0; p_start[4] = p_stop[4] = p_period[4] = 0.0; p_flag[4] = 0; } iarg += 4; } else if (strcmp(arg[iarg],"x") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[0] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = force->numeric(FLERR,arg[iarg+2]); p_period[0] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = 1; deviatoric_flag = 1; iarg += 4; } else if (strcmp(arg[iarg],"y") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[1] = force->numeric(FLERR,arg[iarg+1]); p_stop[1] = force->numeric(FLERR,arg[iarg+2]); p_period[1] = force->numeric(FLERR,arg[iarg+3]); p_flag[1] = 1; deviatoric_flag = 1; iarg += 4; } else if (strcmp(arg[iarg],"z") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[2] = 1; deviatoric_flag = 1; iarg += 4; if (dimension == 2) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); } else if (strcmp(arg[iarg],"yz") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[3] = force->numeric(FLERR,arg[iarg+1]); p_stop[3] = force->numeric(FLERR,arg[iarg+2]); p_period[3] = force->numeric(FLERR,arg[iarg+3]); p_flag[3] = 1; deviatoric_flag = 1; scaleyz = 0; iarg += 4; if (dimension == 2) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); } else if (strcmp(arg[iarg],"xz") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[4] = force->numeric(FLERR,arg[iarg+1]); p_stop[4] = force->numeric(FLERR,arg[iarg+2]); p_period[4] = force->numeric(FLERR,arg[iarg+3]); p_flag[4] = 1; deviatoric_flag = 1; scalexz = 0; iarg += 4; if (dimension == 2) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); } else if (strcmp(arg[iarg],"xy") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[5] = force->numeric(FLERR,arg[iarg+1]); p_stop[5] = force->numeric(FLERR,arg[iarg+2]); p_period[5] = force->numeric(FLERR,arg[iarg+3]); p_flag[5] = 1; deviatoric_flag = 1; scalexy = 0; iarg += 4; } else if (strcmp(arg[iarg],"couple") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ; else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY; else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ; else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ; else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"drag") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); drag = force->numeric(FLERR,arg[iarg+1]); if (drag < 0.0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"dilate") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; else { allremap = 0; delete [] id_dilate; int n = strlen(arg[iarg+1]) + 1; id_dilate = new char[n]; strcpy(id_dilate,arg[iarg+1]); int idilate = group->find(id_dilate); if (idilate == -1) error->all(FLERR,"Fix nvt/npt/nph dilate group ID does not exist"); } iarg += 2; } else if (strcmp(arg[iarg],"tchain") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); mtchain = force->inumeric(FLERR,arg[iarg+1]); // used by FixNVTSllod to preserve non-default value mtchain_default_flag = 0; if (mtchain < 1) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"pchain") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); mpchain = force->inumeric(FLERR,arg[iarg+1]); if (mpchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"mtk") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1; else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"tloop") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); nc_tchain = force->inumeric(FLERR,arg[iarg+1]); if (nc_tchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"ploop") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); nc_pchain = force->inumeric(FLERR,arg[iarg+1]); if (nc_pchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"nreset") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); nreset_h0 = force->inumeric(FLERR,arg[iarg+1]); if (nreset_h0 < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"scalexy") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1; else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"scalexz") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1; else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"scaleyz") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1; else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"flip") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1; else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"fixedpoint") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]); fixedpoint[1] = force->numeric(FLERR,arg[iarg+2]); fixedpoint[2] = force->numeric(FLERR,arg[iarg+3]); iarg += 4; } else error->all(FLERR,"Illegal fix nvt/npt/nph command"); } // error checks if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4])) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); if (dimension == 2 && (scalexz == 1 || scaleyz == 1 )) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); // require periodicity in tensile dimension if (p_flag[0] && domain->xperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension"); if (p_flag[1] && domain->yperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension"); if (p_flag[2] && domain->zperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension"); // require periodicity in 2nd dim of off-diagonal tilt component if (p_flag[3] && domain->zperiodic == 0) error->all(FLERR, "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); if (p_flag[4] && domain->zperiodic == 0) error->all(FLERR, "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); if (p_flag[5] && domain->yperiodic == 0) error->all(FLERR, "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); if (scaleyz == 1 && domain->zperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph " "with yz scaling when z is non-periodic dimension"); if (scalexz == 1 && domain->zperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph " "with xz scaling when z is non-periodic dimension"); if (scalexy == 1 && domain->yperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph " "with xy scaling when y is non-periodic dimension"); if (p_flag[3] && scaleyz == 1) error->all(FLERR,"Cannot use fix nvt/npt/nph with " "both yz dynamics and yz scaling"); if (p_flag[4] && scalexz == 1) error->all(FLERR,"Cannot use fix nvt/npt/nph with " "both xz dynamics and xz scaling"); if (p_flag[5] && scalexy == 1) error->all(FLERR,"Cannot use fix nvt/npt/nph with " "both xy dynamics and xy scaling"); if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5])) error->all(FLERR,"Can not specify Pxy/Pxz/Pyz in " "fix nvt/npt/nph with non-triclinic box"); if (pcouple == XYZ && dimension == 3 && (p_start[0] != p_start[1] || p_start[0] != p_start[2] || p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] || p_period[0] != p_period[1] || p_period[0] != p_period[2])) error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings"); if (pcouple == XYZ && dimension == 2 && (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || p_period[0] != p_period[1])) error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings"); if (pcouple == XY && (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || p_period[0] != p_period[1])) error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings"); if (pcouple == YZ && (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] || p_period[1] != p_period[2])) error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings"); if (pcouple == XZ && (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] || p_period[0] != p_period[2])) error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings"); if ((tstat_flag && t_period <= 0.0) || (p_flag[0] && p_period[0] <= 0.0) || (p_flag[1] && p_period[1] <= 0.0) || (p_flag[2] && p_period[2] <= 0.0) || (p_flag[3] && p_period[3] <= 0.0) || (p_flag[4] && p_period[4] <= 0.0) || (p_flag[5] && p_period[5] <= 0.0)) error->all(FLERR,"Fix nvt/npt/nph damping parameters must be > 0.0"); // set pstat_flag and box change and restart_pbc variables pstat_flag = 0; for (int i = 0; i < 6; i++) if (p_flag[i]) pstat_flag = 1; if (pstat_flag) { if (p_flag[0] || p_flag[1] || p_flag[2]) box_change_size = 1; if (p_flag[3] || p_flag[4] || p_flag[5]) box_change_shape = 1; no_change_box = 1; if (allremap == 0) restart_pbc = 1; } // pstyle = TRICLINIC if any off-diagonal term is controlled -> 6 dof // else pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof // else pstyle = ANISO -> 3 dof if (p_flag[3] || p_flag[4] || p_flag[5]) pstyle = TRICLINIC; else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO; else pstyle = ANISO; // pre_exchange only required if flips can occur due to shape changes pre_exchange_flag = 0; if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5])) pre_exchange_flag = 1; if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0)) pre_exchange_flag = 1; // convert input periods to frequencies t_freq = 0.0; p_freq[0] = p_freq[1] = p_freq[2] = p_freq[3] = p_freq[4] = p_freq[5] = 0.0; if (tstat_flag) t_freq = 1.0 / t_period; if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; if (p_flag[3]) p_freq[3] = 1.0 / p_period[3]; if (p_flag[4]) p_freq[4] = 1.0 / p_period[4]; if (p_flag[5]) p_freq[5] = 1.0 / p_period[5]; // Nose/Hoover temp and pressure init size_vector = 0; if (tstat_flag) { int ich; eta = new double[mtchain]; // add one extra dummy thermostat, set to zero eta_dot = new double[mtchain+1]; eta_dot[mtchain] = 0.0; eta_dotdot = new double[mtchain]; for (ich = 0; ich < mtchain; ich++) { eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0; } eta_mass = new double[mtchain]; size_vector += 2*2*mtchain; } if (pstat_flag) { omega[0] = omega[1] = omega[2] = 0.0; omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0; omega_mass[0] = omega_mass[1] = omega_mass[2] = 0.0; omega[3] = omega[4] = omega[5] = 0.0; omega_dot[3] = omega_dot[4] = omega_dot[5] = 0.0; omega_mass[3] = omega_mass[4] = omega_mass[5] = 0.0; if (pstyle == ISO) size_vector += 2*2*1; else if (pstyle == ANISO) size_vector += 2*2*3; else if (pstyle == TRICLINIC) size_vector += 2*2*6; if (mpchain) { int ich; etap = new double[mpchain]; // add one extra dummy thermostat, set to zero etap_dot = new double[mpchain+1]; etap_dot[mpchain] = 0.0; etap_dotdot = new double[mpchain]; for (ich = 0; ich < mpchain; ich++) { etap[ich] = etap_dot[ich] = etap_dotdot[ich] = 0.0; } etap_mass = new double[mpchain]; size_vector += 2*2*mpchain; } if (deviatoric_flag) size_vector += 1; } nrigid = 0; rfix = NULL; if (pre_exchange_flag) irregular = new Irregular(lmp); else irregular = NULL; // initialize vol0,t0 to zero to signal uninitialized // values then assigned in init(), if necessary vol0 = t0 = 0.0; } /* ---------------------------------------------------------------------- */ FixNH::~FixNH() { + if (copymode) return; + delete [] id_dilate; delete [] rfix; delete irregular; // delete temperature and pressure if fix created them if (tflag) modify->delete_compute(id_temp); delete [] id_temp; if (tstat_flag) { delete [] eta; delete [] eta_dot; delete [] eta_dotdot; delete [] eta_mass; } if (pstat_flag) { if (pflag) modify->delete_compute(id_press); delete [] id_press; if (mpchain) { delete [] etap; delete [] etap_dot; delete [] etap_dotdot; delete [] etap_mass; } } } /* ---------------------------------------------------------------------- */ int FixNH::setmask() { int mask = 0; mask |= INITIAL_INTEGRATE; mask |= FINAL_INTEGRATE; mask |= THERMO_ENERGY; mask |= INITIAL_INTEGRATE_RESPA; mask |= FINAL_INTEGRATE_RESPA; if (pre_exchange_flag) mask |= PRE_EXCHANGE; return mask; } /* ---------------------------------------------------------------------- */ void FixNH::init() { // recheck that dilate group has not been deleted if (allremap == 0) { int idilate = group->find(id_dilate); if (idilate == -1) error->all(FLERR,"Fix nvt/npt/nph dilate group ID does not exist"); dilate_group_bit = group->bitmask[idilate]; } // ensure no conflict with fix deform if (pstat_flag) for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || (p_flag[2] && dimflag[2]) || (p_flag[3] && dimflag[3]) || (p_flag[4] && dimflag[4]) || (p_flag[5] && dimflag[5])) error->all(FLERR,"Cannot use fix npt and fix deform on " "same component of stress tensor"); } // set temperature and pressure ptrs int icompute = modify->find_compute(id_temp); if (icompute < 0) error->all(FLERR,"Temperature ID for fix nvt/npt does not exist"); temperature = modify->compute[icompute]; if (temperature->tempbias) which = BIAS; else which = NOBIAS; if (pstat_flag) { icompute = modify->find_compute(id_press); if (icompute < 0) error->all(FLERR,"Pressure ID for fix npt/nph does not exist"); pressure = modify->compute[icompute]; } // set timesteps and frequencies dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; dthalf = 0.5 * update->dt; dt4 = 0.25 * update->dt; dt8 = 0.125 * update->dt; dto = dthalf; p_freq_max = 0.0; if (pstat_flag) { p_freq_max = MAX(p_freq[0],p_freq[1]); p_freq_max = MAX(p_freq_max,p_freq[2]); if (pstyle == TRICLINIC) { p_freq_max = MAX(p_freq_max,p_freq[3]); p_freq_max = MAX(p_freq_max,p_freq[4]); p_freq_max = MAX(p_freq_max,p_freq[5]); } pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain); } if (tstat_flag) tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain); // tally the number of dimensions that are barostatted // set initial volume and reference cell, if not already done if (pstat_flag) { pdim = p_flag[0] + p_flag[1] + p_flag[2]; if (vol0 == 0.0) { if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd; else vol0 = domain->xprd * domain->yprd; h0_inv[0] = domain->h_inv[0]; h0_inv[1] = domain->h_inv[1]; h0_inv[2] = domain->h_inv[2]; h0_inv[3] = domain->h_inv[3]; h0_inv[4] = domain->h_inv[4]; h0_inv[5] = domain->h_inv[5]; } } boltz = force->boltz; nktv2p = force->nktv2p; if (force->kspace) kspace_flag = 1; else kspace_flag = 0; if (strstr(update->integrate_style,"respa")) { nlevels_respa = ((Respa *) update->integrate)->nlevels; step_respa = ((Respa *) update->integrate)->step; dto = 0.5*step_respa[0]; } // detect if any rigid fixes exist so rigid bodies move when box is remapped // rfix[] = indices to each fix rigid delete [] rfix; nrigid = 0; rfix = NULL; for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) nrigid++; if (nrigid) { rfix = new int[nrigid]; nrigid = 0; for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; } } /* ---------------------------------------------------------------------- compute T,P before integrator starts ------------------------------------------------------------------------- */ void FixNH::setup(int vflag) { // t_target is needed by NPH and NPT in compute_scalar() // If no thermostat or using fix nphug, // t_target must be defined by other means. if (tstat_flag && strcmp(style,"nphug") != 0) { compute_temp_target(); } else if (pstat_flag) { // t0 = reference temperature for masses // cannot be done in init() b/c temperature cannot be called there // is b/c Modify::init() inits computes after fixes due to dof dependence // guesstimate a unit-dependent t0 if actual T = 0.0 // if it was read in from a restart file, leave it be if (t0 == 0.0) { t0 = temperature->compute_scalar(); if (t0 == 0.0) { if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0; else t0 = 300.0; } } t_target = t0; } if (pstat_flag) compute_press_target(); t_current = temperature->compute_scalar(); tdof = temperature->dof; if (pstat_flag) { if (pstyle == ISO) pressure->compute_scalar(); else pressure->compute_vector(); couple(); pressure->addstep(update->ntimestep+1); } // masses and initial forces on thermostat variables if (tstat_flag) { eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); for (int ich = 1; ich < mtchain; ich++) eta_mass[ich] = boltz * t_target / (t_freq*t_freq); for (int ich = 1; ich < mtchain; ich++) { eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] - boltz * t_target) / eta_mass[ich]; } } // masses and initial forces on barostat variables if (pstat_flag) { double kt = boltz * t_target; double nkt = atom->natoms * kt; for (int i = 0; i < 3; i++) if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); if (pstyle == TRICLINIC) { for (int i = 3; i < 6; i++) if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); } // masses and initial forces on barostat thermostat variables if (mpchain) { etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max); for (int ich = 1; ich < mpchain; ich++) etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max); for (int ich = 1; ich < mpchain; ich++) etap_dotdot[ich] = (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - boltz * t_target) / etap_mass[ich]; } } } /* ---------------------------------------------------------------------- 1st half of Verlet update ------------------------------------------------------------------------- */ void FixNH::initial_integrate(int vflag) { // update eta_press_dot if (pstat_flag && mpchain) nhc_press_integrate(); // update eta_dot if (tstat_flag) { compute_temp_target(); nhc_temp_integrate(); } // need to recompute pressure to account for change in KE // t_current is up-to-date, but compute_temperature is not // compute appropriately coupled elements of mvv_current if (pstat_flag) { if (pstyle == ISO) { temperature->compute_scalar(); pressure->compute_scalar(); } else { temperature->compute_vector(); pressure->compute_vector(); } couple(); pressure->addstep(update->ntimestep+1); } if (pstat_flag) { compute_press_target(); nh_omega_dot(); nh_v_press(); } nve_v(); // remap simulation box by 1/2 step if (pstat_flag) remap(); nve_x(); // remap simulation box by 1/2 step // redo KSpace coeffs since volume has changed if (pstat_flag) { remap(); if (kspace_flag) force->kspace->setup(); } } /* ---------------------------------------------------------------------- 2nd half of Verlet update ------------------------------------------------------------------------- */ void FixNH::final_integrate() { nve_v(); // re-compute temp before nh_v_press() // only needed for temperature computes with BIAS on reneighboring steps: // b/c some biases store per-atom values (e.g. temp/profile) // per-atom values are invalid if reneigh/comm occurred // since temp->compute() in initial_integrate() if (which == BIAS && neighbor->ago == 0) t_current = temperature->compute_scalar(); if (pstat_flag) nh_v_press(); // compute new T,P after velocities rescaled by nh_v_press() // compute appropriately coupled elements of mvv_current t_current = temperature->compute_scalar(); tdof = temperature->dof; if (pstat_flag) { if (pstyle == ISO) pressure->compute_scalar(); else pressure->compute_vector(); couple(); pressure->addstep(update->ntimestep+1); } if (pstat_flag) nh_omega_dot(); // update eta_dot // update eta_press_dot if (tstat_flag) nhc_temp_integrate(); if (pstat_flag && mpchain) nhc_press_integrate(); } /* ---------------------------------------------------------------------- */ void FixNH::initial_integrate_respa(int vflag, int ilevel, int iloop) { // set timesteps by level dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; dthalf = 0.5 * step_respa[ilevel]; // outermost level - update eta_dot and omega_dot, apply to v // all other levels - NVE update of v // x,v updates only performed for atoms in group if (ilevel == nlevels_respa-1) { // update eta_press_dot if (pstat_flag && mpchain) nhc_press_integrate(); // update eta_dot if (tstat_flag) { compute_temp_target(); nhc_temp_integrate(); } // recompute pressure to account for change in KE // t_current is up-to-date, but compute_temperature is not // compute appropriately coupled elements of mvv_current if (pstat_flag) { if (pstyle == ISO) { temperature->compute_scalar(); pressure->compute_scalar(); } else { temperature->compute_vector(); pressure->compute_vector(); } couple(); pressure->addstep(update->ntimestep+1); } if (pstat_flag) { compute_press_target(); nh_omega_dot(); nh_v_press(); } nve_v(); } else nve_v(); // innermost level - also update x only for atoms in group // if barostat, perform 1/2 step remap before and after if (ilevel == 0) { if (pstat_flag) remap(); nve_x(); if (pstat_flag) remap(); } // if barostat, redo KSpace coeffs at outermost level, // since volume has changed if (ilevel == nlevels_respa-1 && kspace_flag && pstat_flag) force->kspace->setup(); } /* ---------------------------------------------------------------------- */ void FixNH::final_integrate_respa(int ilevel, int iloop) { // set timesteps by level dtf = 0.5 * step_respa[ilevel] * force->ftm2v; dthalf = 0.5 * step_respa[ilevel]; // outermost level - update eta_dot and omega_dot, apply via final_integrate // all other levels - NVE update of v if (ilevel == nlevels_respa-1) final_integrate(); else nve_v(); } /* ---------------------------------------------------------------------- */ void FixNH::couple() { double *tensor = pressure->vector; if (pstyle == ISO) p_current[0] = p_current[1] = p_current[2] = pressure->scalar; else if (pcouple == XYZ) { double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); p_current[0] = p_current[1] = p_current[2] = ave; } else if (pcouple == XY) { double ave = 0.5 * (tensor[0] + tensor[1]); p_current[0] = p_current[1] = ave; p_current[2] = tensor[2]; } else if (pcouple == YZ) { double ave = 0.5 * (tensor[1] + tensor[2]); p_current[1] = p_current[2] = ave; p_current[0] = tensor[0]; } else if (pcouple == XZ) { double ave = 0.5 * (tensor[0] + tensor[2]); p_current[0] = p_current[2] = ave; p_current[1] = tensor[1]; } else { p_current[0] = tensor[0]; p_current[1] = tensor[1]; p_current[2] = tensor[2]; } // switch order from xy-xz-yz to Voigt if (pstyle == TRICLINIC) { p_current[3] = tensor[5]; p_current[4] = tensor[4]; p_current[5] = tensor[3]; } } /* ---------------------------------------------------------------------- change box size remap all atoms or dilate group atoms depending on allremap flag if rigid bodies exist, scale rigid body centers-of-mass ------------------------------------------------------------------------- */ void FixNH::remap() { int i; double oldlo,oldhi; double expfac; double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; double *h = domain->h; // omega is not used, except for book-keeping for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i]; // convert pertinent atoms and rigid bodies to lamda coords if (allremap) domain->x2lamda(nlocal); else { for (i = 0; i < nlocal; i++) if (mask[i] & dilate_group_bit) domain->x2lamda(x[i],x[i]); } if (nrigid) for (i = 0; i < nrigid; i++) modify->fix[rfix[i]]->deform(0); // reset global and local box to new size/shape // this operation corresponds to applying the // translate and scale operations // corresponding to the solution of the following ODE: // // h_dot = omega_dot * h // // where h_dot, omega_dot and h are all upper-triangular // 3x3 tensors. In Voigt notation, the elements of the // RHS product tensor are: // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1] // // Ordering of operations preserves time symmetry. double dto2 = dto/2.0; double dto4 = dto/4.0; double dto8 = dto/8.0; // off-diagonal components, first half if (pstyle == TRICLINIC) { if (p_flag[4]) { expfac = exp(dto8*omega_dot[0]); h[4] *= expfac; h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); h[4] *= expfac; } if (p_flag[3]) { expfac = exp(dto4*omega_dot[1]); h[3] *= expfac; h[3] += dto2*(omega_dot[3]*h[2]); h[3] *= expfac; } if (p_flag[5]) { expfac = exp(dto4*omega_dot[0]); h[5] *= expfac; h[5] += dto2*(omega_dot[5]*h[1]); h[5] *= expfac; } if (p_flag[4]) { expfac = exp(dto8*omega_dot[0]); h[4] *= expfac; h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); h[4] *= expfac; } } // scale diagonal components // scale tilt factors with cell, if set if (p_flag[0]) { oldlo = domain->boxlo[0]; oldhi = domain->boxhi[0]; expfac = exp(dto*omega_dot[0]); domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0]; domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0]; } if (p_flag[1]) { oldlo = domain->boxlo[1]; oldhi = domain->boxhi[1]; expfac = exp(dto*omega_dot[1]); domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1]; domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1]; if (scalexy) h[5] *= expfac; } if (p_flag[2]) { oldlo = domain->boxlo[2]; oldhi = domain->boxhi[2]; expfac = exp(dto*omega_dot[2]); domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2]; domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2]; if (scalexz) h[4] *= expfac; if (scaleyz) h[3] *= expfac; } // off-diagonal components, second half if (pstyle == TRICLINIC) { if (p_flag[4]) { expfac = exp(dto8*omega_dot[0]); h[4] *= expfac; h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); h[4] *= expfac; } if (p_flag[3]) { expfac = exp(dto4*omega_dot[1]); h[3] *= expfac; h[3] += dto2*(omega_dot[3]*h[2]); h[3] *= expfac; } if (p_flag[5]) { expfac = exp(dto4*omega_dot[0]); h[5] *= expfac; h[5] += dto2*(omega_dot[5]*h[1]); h[5] *= expfac; } if (p_flag[4]) { expfac = exp(dto8*omega_dot[0]); h[4] *= expfac; h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); h[4] *= expfac; } } domain->yz = h[3]; domain->xz = h[4]; domain->xy = h[5]; // tilt factor to cell length ratio can not exceed TILTMAX in one step if (domain->yz < -TILTMAX*domain->yprd || domain->yz > TILTMAX*domain->yprd || domain->xz < -TILTMAX*domain->xprd || domain->xz > TILTMAX*domain->xprd || domain->xy < -TILTMAX*domain->xprd || domain->xy > TILTMAX*domain->xprd) error->all(FLERR,"Fix npt/nph has tilted box too far in one step - " "periodic cell is too far from equilibrium state"); domain->set_global_box(); domain->set_local_box(); // convert pertinent atoms and rigid bodies back to box coords if (allremap) domain->lamda2x(nlocal); else { for (i = 0; i < nlocal; i++) if (mask[i] & dilate_group_bit) domain->lamda2x(x[i],x[i]); } if (nrigid) for (i = 0; i < nrigid; i++) modify->fix[rfix[i]]->deform(1); } /* ---------------------------------------------------------------------- pack entire state of Fix into one write ------------------------------------------------------------------------- */ void FixNH::write_restart(FILE *fp) { int nsize = size_restart_global(); double *list; memory->create(list,nsize,"nh:list"); pack_restart_data(list); if (comm->me == 0) { int size = nsize * sizeof(double); fwrite(&size,sizeof(int),1,fp); fwrite(list,sizeof(double),nsize,fp); } memory->destroy(list); } /* ---------------------------------------------------------------------- calculate the number of data to be packed ------------------------------------------------------------------------- */ int FixNH::size_restart_global() { int nsize = 2; if (tstat_flag) nsize += 1 + 2*mtchain; if (pstat_flag) { nsize += 16 + 2*mpchain; if (deviatoric_flag) nsize += 6; } return nsize; } /* ---------------------------------------------------------------------- pack restart data ------------------------------------------------------------------------- */ int FixNH::pack_restart_data(double *list) { int n = 0; list[n++] = tstat_flag; if (tstat_flag) { list[n++] = mtchain; for (int ich = 0; ich < mtchain; ich++) list[n++] = eta[ich]; for (int ich = 0; ich < mtchain; ich++) list[n++] = eta_dot[ich]; } list[n++] = pstat_flag; if (pstat_flag) { list[n++] = omega[0]; list[n++] = omega[1]; list[n++] = omega[2]; list[n++] = omega[3]; list[n++] = omega[4]; list[n++] = omega[5]; list[n++] = omega_dot[0]; list[n++] = omega_dot[1]; list[n++] = omega_dot[2]; list[n++] = omega_dot[3]; list[n++] = omega_dot[4]; list[n++] = omega_dot[5]; list[n++] = vol0; list[n++] = t0; list[n++] = mpchain; if (mpchain) { for (int ich = 0; ich < mpchain; ich++) list[n++] = etap[ich]; for (int ich = 0; ich < mpchain; ich++) list[n++] = etap_dot[ich]; } list[n++] = deviatoric_flag; if (deviatoric_flag) { list[n++] = h0_inv[0]; list[n++] = h0_inv[1]; list[n++] = h0_inv[2]; list[n++] = h0_inv[3]; list[n++] = h0_inv[4]; list[n++] = h0_inv[5]; } } return n; } /* ---------------------------------------------------------------------- use state info from restart file to restart the Fix ------------------------------------------------------------------------- */ void FixNH::restart(char *buf) { int n = 0; double *list = (double *) buf; int flag = static_cast (list[n++]); if (flag) { int m = static_cast (list[n++]); if (tstat_flag && m == mtchain) { for (int ich = 0; ich < mtchain; ich++) eta[ich] = list[n++]; for (int ich = 0; ich < mtchain; ich++) eta_dot[ich] = list[n++]; } else n += 2*m; } flag = static_cast (list[n++]); if (flag) { omega[0] = list[n++]; omega[1] = list[n++]; omega[2] = list[n++]; omega[3] = list[n++]; omega[4] = list[n++]; omega[5] = list[n++]; omega_dot[0] = list[n++]; omega_dot[1] = list[n++]; omega_dot[2] = list[n++]; omega_dot[3] = list[n++]; omega_dot[4] = list[n++]; omega_dot[5] = list[n++]; vol0 = list[n++]; t0 = list[n++]; int m = static_cast (list[n++]); if (pstat_flag && m == mpchain) { for (int ich = 0; ich < mpchain; ich++) etap[ich] = list[n++]; for (int ich = 0; ich < mpchain; ich++) etap_dot[ich] = list[n++]; } else n+=2*m; flag = static_cast (list[n++]); if (flag) { h0_inv[0] = list[n++]; h0_inv[1] = list[n++]; h0_inv[2] = list[n++]; h0_inv[3] = list[n++]; h0_inv[4] = list[n++]; h0_inv[5] = list[n++]; } } } /* ---------------------------------------------------------------------- */ int FixNH::modify_param(int narg, char **arg) { if (strcmp(arg[0],"temp") == 0) { if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); if (tflag) { modify->delete_compute(id_temp); tflag = 0; } delete [] id_temp; int n = strlen(arg[1]) + 1; id_temp = new char[n]; strcpy(id_temp,arg[1]); int icompute = modify->find_compute(arg[1]); if (icompute < 0) error->all(FLERR,"Could not find fix_modify temperature ID"); temperature = modify->compute[icompute]; if (temperature->tempflag == 0) error->all(FLERR, "Fix_modify temperature ID does not compute temperature"); if (temperature->igroup != 0 && comm->me == 0) error->warning(FLERR,"Temperature for fix modify is not for group all"); // reset id_temp of pressure to new temperature ID if (pstat_flag) { icompute = modify->find_compute(id_press); if (icompute < 0) error->all(FLERR,"Pressure ID for fix modify does not exist"); modify->compute[icompute]->reset_extra_compute_fix(id_temp); } return 2; } else if (strcmp(arg[0],"press") == 0) { if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command"); if (pflag) { modify->delete_compute(id_press); pflag = 0; } delete [] id_press; int n = strlen(arg[1]) + 1; id_press = new char[n]; strcpy(id_press,arg[1]); int icompute = modify->find_compute(arg[1]); if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); pressure = modify->compute[icompute]; if (pressure->pressflag == 0) error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); return 2; } return 0; } /* ---------------------------------------------------------------------- */ double FixNH::compute_scalar() { int i; double volume; double energy; double kt = boltz * t_target; double lkt_press = kt; int ich; if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; else volume = domain->xprd * domain->yprd; energy = 0.0; // thermostat chain energy is equivalent to Eq. (2) in // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117 // Sum(0.5*p_eta_k^2/Q_k,k=1,M) + L*k*T*eta_1 + Sum(k*T*eta_k,k=2,M), // where L = tdof // M = mtchain // p_eta_k = Q_k*eta_dot[k-1] // Q_1 = L*k*T/t_freq^2 // Q_k = k*T/t_freq^2, k > 1 if (tstat_flag) { energy += ke_target * eta[0] + 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0]; for (ich = 1; ich < mtchain; ich++) energy += kt * eta[ich] + 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich]; } // barostat energy is equivalent to Eq. (8) in // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117 // Sum(0.5*p_omega^2/W + P*V), // where N = natoms // p_omega = W*omega_dot // W = N*k*T/p_freq^2 // sum is over barostatted dimensions if (pstat_flag) { for (i = 0; i < 3; i++) if (p_flag[i]) energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] + p_hydro*(volume-vol0) / (pdim*nktv2p); if (pstyle == TRICLINIC) { for (i = 3; i < 6; i++) if (p_flag[i]) energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i]; } // extra contributions from thermostat chain for barostat if (mpchain) { energy += lkt_press * etap[0] + 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0]; for (ich = 1; ich < mpchain; ich++) energy += kt * etap[ich] + 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich]; } // extra contribution from strain energy if (deviatoric_flag) energy += compute_strain_energy(); } return energy; } /* ---------------------------------------------------------------------- return a single element of the following vectors, in this order: eta[tchain], eta_dot[tchain], omega[ndof], omega_dot[ndof] etap[pchain], etap_dot[pchain], PE_eta[tchain], KE_eta_dot[tchain] PE_omega[ndof], KE_omega_dot[ndof], PE_etap[pchain], KE_etap_dot[pchain] PE_strain[1] if no thermostat exists, related quantities are omitted from the list if no barostat exists, related quantities are omitted from the list ndof = 1,3,6 degrees of freedom for pstyle = ISO,ANISO,TRI ------------------------------------------------------------------------- */ double FixNH::compute_vector(int n) { int ilen; if (tstat_flag) { ilen = mtchain; if (n < ilen) return eta[n]; n -= ilen; ilen = mtchain; if (n < ilen) return eta_dot[n]; n -= ilen; } if (pstat_flag) { if (pstyle == ISO) { ilen = 1; if (n < ilen) return omega[n]; n -= ilen; } else if (pstyle == ANISO) { ilen = 3; if (n < ilen) return omega[n]; n -= ilen; } else { ilen = 6; if (n < ilen) return omega[n]; n -= ilen; } if (pstyle == ISO) { ilen = 1; if (n < ilen) return omega_dot[n]; n -= ilen; } else if (pstyle == ANISO) { ilen = 3; if (n < ilen) return omega_dot[n]; n -= ilen; } else { ilen = 6; if (n < ilen) return omega_dot[n]; n -= ilen; } if (mpchain) { ilen = mpchain; if (n < ilen) return etap[n]; n -= ilen; ilen = mpchain; if (n < ilen) return etap_dot[n]; n -= ilen; } } double volume; double kt = boltz * t_target; double lkt_press = kt; int ich; if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; else volume = domain->xprd * domain->yprd; if (tstat_flag) { ilen = mtchain; if (n < ilen) { ich = n; if (ich == 0) return ke_target * eta[0]; else return kt * eta[ich]; } n -= ilen; ilen = mtchain; if (n < ilen) { ich = n; if (ich == 0) return 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0]; else return 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich]; } n -= ilen; } if (pstat_flag) { if (pstyle == ISO) { ilen = 1; if (n < ilen) return p_hydro*(volume-vol0) / nktv2p; n -= ilen; } else if (pstyle == ANISO) { ilen = 3; if (n < ilen) { if (p_flag[n]) return p_hydro*(volume-vol0) / (pdim*nktv2p); else return 0.0; } n -= ilen; } else { ilen = 6; if (n < ilen) { if (n > 2) return 0.0; else if (p_flag[n]) return p_hydro*(volume-vol0) / (pdim*nktv2p); else return 0.0; } n -= ilen; } if (pstyle == ISO) { ilen = 1; if (n < ilen) return pdim*0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; n -= ilen; } else if (pstyle == ANISO) { ilen = 3; if (n < ilen) { if (p_flag[n]) return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; else return 0.0; } n -= ilen; } else { ilen = 6; if (n < ilen) { if (p_flag[n]) return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; else return 0.0; } n -= ilen; } if (mpchain) { ilen = mpchain; if (n < ilen) { ich = n; if (ich == 0) return lkt_press * etap[0]; else return kt * etap[ich]; } n -= ilen; ilen = mpchain; if (n < ilen) { ich = n; if (ich == 0) return 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0]; else return 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich]; } n -= ilen; } if (deviatoric_flag) { ilen = 1; if (n < ilen) return compute_strain_energy(); n -= ilen; } } return 0.0; } /* ---------------------------------------------------------------------- */ void FixNH::reset_target(double t_new) { t_target = t_start = t_stop = t_new; } /* ---------------------------------------------------------------------- */ void FixNH::reset_dt() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; dthalf = 0.5 * update->dt; dt4 = 0.25 * update->dt; dt8 = 0.125 * update->dt; dto = dthalf; // If using respa, then remap is performed in innermost level if (strstr(update->integrate_style,"respa")) dto = 0.5*step_respa[0]; if (pstat_flag) pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain); if (tstat_flag) tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain); } /* ---------------------------------------------------------------------- extract thermostat properties ------------------------------------------------------------------------- */ void *FixNH::extract(const char *str, int &dim) { dim=0; if (strcmp(str,"t_target") == 0) { return &t_target; } else if (strcmp(str,"mtchain") == 0) { return &mtchain; } dim=1; if (strcmp(str,"eta") == 0) { return η } return NULL; } /* ---------------------------------------------------------------------- perform half-step update of chain thermostat variables ------------------------------------------------------------------------- */ void FixNH::nhc_temp_integrate() { int ich; double expfac; double kecurrent = tdof * boltz * t_current; // Update masses, to preserve initial freq, if flag set if (eta_mass_flag) { eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); for (int ich = 1; ich < mtchain; ich++) eta_mass[ich] = boltz * t_target / (t_freq*t_freq); } if (eta_mass[0] > 0.0) eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0]; else eta_dotdot[0] = 0.0; double ncfac = 1.0/nc_tchain; for (int iloop = 0; iloop < nc_tchain; iloop++) { for (ich = mtchain-1; ich > 0; ich--) { expfac = exp(-ncfac*dt8*eta_dot[ich+1]); eta_dot[ich] *= expfac; eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4; eta_dot[ich] *= tdrag_factor; eta_dot[ich] *= expfac; } expfac = exp(-ncfac*dt8*eta_dot[1]); eta_dot[0] *= expfac; eta_dot[0] += eta_dotdot[0] * ncfac*dt4; eta_dot[0] *= tdrag_factor; eta_dot[0] *= expfac; factor_eta = exp(-ncfac*dthalf*eta_dot[0]); nh_v_temp(); // rescale temperature due to velocity scaling // should not be necessary to explicitly recompute the temperature t_current *= factor_eta*factor_eta; kecurrent = tdof * boltz * t_current; if (eta_mass[0] > 0.0) eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0]; else eta_dotdot[0] = 0.0; for (ich = 0; ich < mtchain; ich++) eta[ich] += ncfac*dthalf*eta_dot[ich]; eta_dot[0] *= expfac; eta_dot[0] += eta_dotdot[0] * ncfac*dt4; eta_dot[0] *= expfac; for (ich = 1; ich < mtchain; ich++) { expfac = exp(-ncfac*dt8*eta_dot[ich+1]); eta_dot[ich] *= expfac; eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] - boltz * t_target)/eta_mass[ich]; eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4; eta_dot[ich] *= expfac; } } } /* ---------------------------------------------------------------------- perform half-step update of chain thermostat variables for barostat scale barostat velocities ------------------------------------------------------------------------- */ void FixNH::nhc_press_integrate() { int ich,i; double expfac,factor_etap,kecurrent; double kt = boltz * t_target; double lkt_press = kt; // Update masses, to preserve initial freq, if flag set if (omega_mass_flag) { double nkt = atom->natoms * kt; for (int i = 0; i < 3; i++) if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); if (pstyle == TRICLINIC) { for (int i = 3; i < 6; i++) if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); } } if (etap_mass_flag) { if (mpchain) { etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max); for (int ich = 1; ich < mpchain; ich++) etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max); for (int ich = 1; ich < mpchain; ich++) etap_dotdot[ich] = (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - boltz * t_target) / etap_mass[ich]; } } kecurrent = 0.0; for (i = 0; i < 3; i++) if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; if (pstyle == TRICLINIC) { for (i = 3; i < 6; i++) if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; } etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0]; double ncfac = 1.0/nc_pchain; for (int iloop = 0; iloop < nc_pchain; iloop++) { for (ich = mpchain-1; ich > 0; ich--) { expfac = exp(-ncfac*dt8*etap_dot[ich+1]); etap_dot[ich] *= expfac; etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4; etap_dot[ich] *= pdrag_factor; etap_dot[ich] *= expfac; } expfac = exp(-ncfac*dt8*etap_dot[1]); etap_dot[0] *= expfac; etap_dot[0] += etap_dotdot[0] * ncfac*dt4; etap_dot[0] *= pdrag_factor; etap_dot[0] *= expfac; for (ich = 0; ich < mpchain; ich++) etap[ich] += ncfac*dthalf*etap_dot[ich]; factor_etap = exp(-ncfac*dthalf*etap_dot[0]); for (i = 0; i < 3; i++) if (p_flag[i]) omega_dot[i] *= factor_etap; if (pstyle == TRICLINIC) { for (i = 3; i < 6; i++) if (p_flag[i]) omega_dot[i] *= factor_etap; } kecurrent = 0.0; for (i = 0; i < 3; i++) if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; if (pstyle == TRICLINIC) { for (i = 3; i < 6; i++) if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; } etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0]; etap_dot[0] *= expfac; etap_dot[0] += etap_dotdot[0] * ncfac*dt4; etap_dot[0] *= expfac; for (ich = 1; ich < mpchain; ich++) { expfac = exp(-ncfac*dt8*etap_dot[ich+1]); etap_dot[ich] *= expfac; etap_dotdot[ich] = (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - boltz*t_target) / etap_mass[ich]; etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4; etap_dot[ich] *= expfac; } } } /* ---------------------------------------------------------------------- perform half-step barostat scaling of velocities -----------------------------------------------------------------------*/ void FixNH::nh_v_press() { double factor[3]; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; factor[0] = exp(-dt4*(omega_dot[0]+mtk_term2)); factor[1] = exp(-dt4*(omega_dot[1]+mtk_term2)); factor[2] = exp(-dt4*(omega_dot[2]+mtk_term2)); if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { v[i][0] *= factor[0]; v[i][1] *= factor[1]; v[i][2] *= factor[2]; if (pstyle == TRICLINIC) { v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]); v[i][1] += -dthalf*v[i][2]*omega_dot[3]; } v[i][0] *= factor[0]; v[i][1] *= factor[1]; v[i][2] *= factor[2]; } } } else if (which == BIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); v[i][0] *= factor[0]; v[i][1] *= factor[1]; v[i][2] *= factor[2]; if (pstyle == TRICLINIC) { v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]); v[i][1] += -dthalf*v[i][2]*omega_dot[3]; } v[i][0] *= factor[0]; v[i][1] *= factor[1]; v[i][2] *= factor[2]; temperature->restore_bias(i,v[i]); } } } } /* ---------------------------------------------------------------------- perform half-step update of velocities -----------------------------------------------------------------------*/ void FixNH::nve_v() { double dtfm; double **v = atom->v; double **f = atom->f; double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; if (rmass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { dtfm = dtf / rmass[i]; v[i][0] += dtfm*f[i][0]; v[i][1] += dtfm*f[i][1]; v[i][2] += dtfm*f[i][2]; } } } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { dtfm = dtf / mass[type[i]]; v[i][0] += dtfm*f[i][0]; v[i][1] += dtfm*f[i][1]; v[i][2] += dtfm*f[i][2]; } } } } /* ---------------------------------------------------------------------- perform full-step update of positions -----------------------------------------------------------------------*/ void FixNH::nve_x() { double **x = atom->x; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; // x update by full step only for atoms in group for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { x[i][0] += dtv * v[i][0]; x[i][1] += dtv * v[i][1]; x[i][2] += dtv * v[i][2]; } } } /* ---------------------------------------------------------------------- perform half-step thermostat scaling of velocities -----------------------------------------------------------------------*/ void FixNH::nh_v_temp() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { v[i][0] *= factor_eta; v[i][1] *= factor_eta; v[i][2] *= factor_eta; } } } else if (which == BIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); v[i][0] *= factor_eta; v[i][1] *= factor_eta; v[i][2] *= factor_eta; temperature->restore_bias(i,v[i]); } } } } /* ---------------------------------------------------------------------- compute sigma tensor needed whenever p_target or h0_inv changes -----------------------------------------------------------------------*/ void FixNH::compute_sigma() { // if nreset_h0 > 0, reset vol0 and h0_inv // every nreset_h0 timesteps if (nreset_h0 > 0) { int delta = update->ntimestep - update->beginstep; if (delta % nreset_h0 == 0) { if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd; else vol0 = domain->xprd * domain->yprd; h0_inv[0] = domain->h_inv[0]; h0_inv[1] = domain->h_inv[1]; h0_inv[2] = domain->h_inv[2]; h0_inv[3] = domain->h_inv[3]; h0_inv[4] = domain->h_inv[4]; h0_inv[5] = domain->h_inv[5]; } } // generate upper-triangular half of // sigma = vol0*h0inv*(p_target-p_hydro)*h0inv^t // units of sigma are are PV/L^2 e.g. atm.A // // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ] // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ] // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ] sigma[0] = vol0*(h0_inv[0]*((p_target[0]-p_hydro)*h0_inv[0] + p_target[5]*h0_inv[5]+p_target[4]*h0_inv[4]) + h0_inv[5]*(p_target[5]*h0_inv[0] + (p_target[1]-p_hydro)*h0_inv[5]+p_target[3]*h0_inv[4]) + h0_inv[4]*(p_target[4]*h0_inv[0]+p_target[3]*h0_inv[5] + (p_target[2]-p_hydro)*h0_inv[4])); sigma[1] = vol0*(h0_inv[1]*((p_target[1]-p_hydro)*h0_inv[1] + p_target[3]*h0_inv[3]) + h0_inv[3]*(p_target[3]*h0_inv[1] + (p_target[2]-p_hydro)*h0_inv[3])); sigma[2] = vol0*(h0_inv[2]*((p_target[2]-p_hydro)*h0_inv[2])); sigma[3] = vol0*(h0_inv[1]*(p_target[3]*h0_inv[2]) + h0_inv[3]*((p_target[2]-p_hydro)*h0_inv[2])); sigma[4] = vol0*(h0_inv[0]*(p_target[4]*h0_inv[2]) + h0_inv[5]*(p_target[3]*h0_inv[2]) + h0_inv[4]*((p_target[2]-p_hydro)*h0_inv[2])); sigma[5] = vol0*(h0_inv[0]*(p_target[5]*h0_inv[1]+p_target[4]*h0_inv[3]) + h0_inv[5]*((p_target[1]-p_hydro)*h0_inv[1]+p_target[3]*h0_inv[3]) + h0_inv[4]*(p_target[3]*h0_inv[1]+(p_target[2]-p_hydro)*h0_inv[3])); } /* ---------------------------------------------------------------------- compute strain energy -----------------------------------------------------------------------*/ double FixNH::compute_strain_energy() { // compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units double* h = domain->h; double d0,d1,d2; d0 = sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) + sigma[5]*( h[1]*h[5]+h[3]*h[4]) + sigma[4]*( h[2]*h[4]); d1 = sigma[5]*( h[5]*h[1]+h[4]*h[3]) + sigma[1]*( h[1]*h[1]+h[3]*h[3]) + sigma[3]*( h[2]*h[3]); d2 = sigma[4]*( h[4]*h[2]) + sigma[3]*( h[3]*h[2]) + sigma[2]*( h[2]*h[2]); double energy = 0.5*(d0+d1+d2)/nktv2p; return energy; } /* ---------------------------------------------------------------------- compute deviatoric barostat force = h*sigma*h^t -----------------------------------------------------------------------*/ void FixNH::compute_deviatoric() { // generate upper-triangular part of h*sigma*h^t // units of fdev are are PV, e.g. atm*A^3 // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ] // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ] // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ] double* h = domain->h; fdev[0] = h[0]*(sigma[0]*h[0]+sigma[5]*h[5]+sigma[4]*h[4]) + h[5]*(sigma[5]*h[0]+sigma[1]*h[5]+sigma[3]*h[4]) + h[4]*(sigma[4]*h[0]+sigma[3]*h[5]+sigma[2]*h[4]); fdev[1] = h[1]*( sigma[1]*h[1]+sigma[3]*h[3]) + h[3]*( sigma[3]*h[1]+sigma[2]*h[3]); fdev[2] = h[2]*( sigma[2]*h[2]); fdev[3] = h[1]*( sigma[3]*h[2]) + h[3]*( sigma[2]*h[2]); fdev[4] = h[0]*( sigma[4]*h[2]) + h[5]*( sigma[3]*h[2]) + h[4]*( sigma[2]*h[2]); fdev[5] = h[0]*( sigma[5]*h[1]+sigma[4]*h[3]) + h[5]*( sigma[1]*h[1]+sigma[3]*h[3]) + h[4]*( sigma[3]*h[1]+sigma[2]*h[3]); } /* ---------------------------------------------------------------------- compute target temperature and kinetic energy -----------------------------------------------------------------------*/ void FixNH::compute_temp_target() { double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; t_target = t_start + delta * (t_stop-t_start); ke_target = tdof * boltz * t_target; } /* ---------------------------------------------------------------------- compute hydrostatic target pressure -----------------------------------------------------------------------*/ void FixNH::compute_press_target() { double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; p_hydro = 0.0; for (int i = 0; i < 3; i++) if (p_flag[i]) { p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); p_hydro += p_target[i]; } p_hydro /= pdim; if (pstyle == TRICLINIC) for (int i = 3; i < 6; i++) p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); // if deviatoric, recompute sigma each time p_target changes if (deviatoric_flag) compute_sigma(); } /* ---------------------------------------------------------------------- update omega_dot, omega -----------------------------------------------------------------------*/ void FixNH::nh_omega_dot() { double f_omega,volume; if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; else volume = domain->xprd*domain->yprd; if (deviatoric_flag) compute_deviatoric(); mtk_term1 = 0.0; if (mtk_flag) { if (pstyle == ISO) { mtk_term1 = tdof * boltz * t_current; mtk_term1 /= pdim * atom->natoms; } else { double *mvv_current = temperature->vector; for (int i = 0; i < 3; i++) if (p_flag[i]) mtk_term1 += mvv_current[i]; mtk_term1 /= pdim * atom->natoms; } } for (int i = 0; i < 3; i++) if (p_flag[i]) { f_omega = (p_current[i]-p_hydro)*volume / (omega_mass[i] * nktv2p) + mtk_term1 / omega_mass[i]; if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p); omega_dot[i] += f_omega*dthalf; omega_dot[i] *= pdrag_factor; } mtk_term2 = 0.0; if (mtk_flag) { for (int i = 0; i < 3; i++) if (p_flag[i]) mtk_term2 += omega_dot[i]; mtk_term2 /= pdim * atom->natoms; } if (pstyle == TRICLINIC) { for (int i = 3; i < 6; i++) { if (p_flag[i]) { f_omega = p_current[i]*volume/(omega_mass[i] * nktv2p); if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p); omega_dot[i] += f_omega*dthalf; omega_dot[i] *= pdrag_factor; } } } } /* ---------------------------------------------------------------------- if any tilt ratios exceed limits, set flip = 1 and compute new tilt values do not flip in x or y if non-periodic (can tilt but not flip) this is b/c the box length would be changed (dramatically) by flip if yz tilt exceeded, adjust C vector by one B vector if xz tilt exceeded, adjust C vector by one A vector if xy tilt exceeded, adjust B vector by one A vector check yz first since it may change xz, then xz check comes after if any flip occurs, create new box in domain image_flip() adjusts image flags due to box shape change induced by flip remap() puts atoms outside the new box back into the new box perform irregular on atoms in lamda coords to migrate atoms to new procs important that image_flip comes before remap, since remap may change image flags to new values, making eqs in doc of Domain:image_flip incorrect ------------------------------------------------------------------------- */ void FixNH::pre_exchange() { double xprd = domain->xprd; double yprd = domain->yprd; // flip is only triggered when tilt exceeds 0.5 by DELTAFLIP // this avoids immediate re-flipping due to tilt oscillations double xtiltmax = (0.5+DELTAFLIP)*xprd; double ytiltmax = (0.5+DELTAFLIP)*yprd; int flipxy,flipxz,flipyz; flipxy = flipxz = flipyz = 0; if (domain->yperiodic) { if (domain->yz < -ytiltmax) { domain->yz += yprd; domain->xz += domain->xy; flipyz = 1; } else if (domain->yz >= ytiltmax) { domain->yz -= yprd; domain->xz -= domain->xy; flipyz = -1; } } if (domain->xperiodic) { if (domain->xz < -xtiltmax) { domain->xz += xprd; flipxz = 1; } else if (domain->xz >= xtiltmax) { domain->xz -= xprd; flipxz = -1; } if (domain->xy < -xtiltmax) { domain->xy += xprd; flipxy = 1; } else if (domain->xy >= xtiltmax) { domain->xy -= xprd; flipxy = -1; } } int flip = 0; if (flipxy || flipxz || flipyz) flip = 1; if (flip) { domain->set_global_box(); domain->set_local_box(); domain->image_flip(flipxy,flipxz,flipyz); double **x = atom->x; imageint *image = atom->image; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); domain->x2lamda(atom->nlocal); irregular->migrate_atoms(); domain->lamda2x(atom->nlocal); } } /* ---------------------------------------------------------------------- memory usage of Irregular ------------------------------------------------------------------------- */ double FixNH::memory_usage() { double bytes = 0.0; if (irregular) bytes += irregular->memory_usage(); return bytes; } diff --git a/src/fix_nh.h b/src/fix_nh.h index 71ea86eeb..be645d0be 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -1,260 +1,260 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifndef LMP_FIX_NH_H #define LMP_FIX_NH_H #include "fix.h" namespace LAMMPS_NS { class FixNH : public Fix { public: FixNH(class LAMMPS *, int, char **); virtual ~FixNH(); int setmask(); virtual void init(); virtual void setup(int); virtual void initial_integrate(int); virtual void final_integrate(); void initial_integrate_respa(int, int, int); void final_integrate_respa(int, int); - void pre_exchange(); + virtual void pre_exchange(); double compute_scalar(); virtual double compute_vector(int); void write_restart(FILE *); virtual int pack_restart_data(double *); // pack restart data virtual void restart(char *); int modify_param(int, char **); void reset_target(double); void reset_dt(); virtual void *extract(const char*,int &); double memory_usage(); protected: int dimension,which; double dtv,dtf,dthalf,dt4,dt8,dto; double boltz,nktv2p,tdof; double vol0; // reference volume double t0; // reference temperature // used for barostat mass double t_start,t_stop; double t_current,t_target,ke_target; double t_freq; int tstat_flag; // 1 if control T int pstat_flag; // 1 if control P int pstyle,pcouple,allremap; int p_flag[6]; // 1 if control P on this dim, 0 if not double p_start[6],p_stop[6]; double p_freq[6],p_target[6]; double omega[6],omega_dot[6]; double omega_mass[6]; double p_current[6]; double drag,tdrag_factor; // drag factor on particle thermostat double pdrag_factor; // drag factor on barostat int kspace_flag; // 1 if KSpace invoked, 0 if not int nrigid; // number of rigid fixes int dilate_group_bit; // mask for dilation group int *rfix; // indices of rigid fixes char *id_dilate; // group name to dilate class Irregular *irregular; // for migrating atoms after box flips int nlevels_respa; double *step_respa; char *id_temp,*id_press; class Compute *temperature,*pressure; int tflag,pflag; double *eta,*eta_dot; // chain thermostat for particles double *eta_dotdot; double *eta_mass; int mtchain; // length of chain int mtchain_default_flag; // 1 = mtchain is default double *etap; // chain thermostat for barostat double *etap_dot; double *etap_dotdot; double *etap_mass; int mpchain; // length of chain int mtk_flag; // 0 if using Hoover barostat int pdim; // number of barostatted dims double p_freq_max; // maximum barostat frequency double p_hydro; // hydrostatic target pressure int nc_tchain,nc_pchain; double factor_eta; double sigma[6]; // scaled target stress double fdev[6]; // deviatoric force on barostat int deviatoric_flag; // 0 if target stress tensor is hydrostatic double h0_inv[6]; // h_inv of reference (zero strain) box int nreset_h0; // interval for resetting h0 double mtk_term1,mtk_term2; // Martyna-Tobias-Klein corrections int eta_mass_flag; // 1 if eta_mass updated, 0 if not. int omega_mass_flag; // 1 if omega_mass updated, 0 if not. int etap_mass_flag; // 1 if etap_mass updated, 0 if not. int scaleyz; // 1 if yz scaled with lz int scalexz; // 1 if xz scaled with lz int scalexy; // 1 if xy scaled with ly int flipflag; // 1 if box flips are invoked as needed int pre_exchange_flag; // set if pre_exchange needed for box flips double fixedpoint[3]; // location of dilation fixed-point void couple(); - void remap(); + virtual void remap(); void nhc_temp_integrate(); void nhc_press_integrate(); virtual void nve_x(); // may be overwritten by child classes virtual void nve_v(); virtual void nh_v_press(); virtual void nh_v_temp(); virtual void compute_temp_target(); virtual int size_restart_global(); void compute_sigma(); void compute_deviatoric(); double compute_strain_energy(); void compute_press_target(); void nh_omega_dot(); }; } #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Target temperature for fix nvt/npt/nph cannot be 0.0 Self-explanatory. E: Invalid fix nvt/npt/nph command for a 2d simulation Cannot control z dimension in a 2d model. E: Fix nvt/npt/nph dilate group ID does not exist Self-explanatory. E: Invalid fix nvt/npt/nph command pressure settings If multiple dimensions are coupled, those dimensions must be specified. E: Cannot use fix nvt/npt/nph on a non-periodic dimension When specifying a diagonal pressure component, the dimension must be periodic. E: Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension When specifying an off-diagonal pressure component, the 2nd of the two dimensions must be periodic. E.g. if the xy component is specified, then the y dimension must be periodic. E: Cannot use fix nvt/npt/nph with yz scaling when z is non-periodic dimension The 2nd dimension in the barostatted tilt factor must be periodic. E: Cannot use fix nvt/npt/nph with xz scaling when z is non-periodic dimension The 2nd dimension in the barostatted tilt factor must be periodic. E: Cannot use fix nvt/npt/nph with xy scaling when y is non-periodic dimension The 2nd dimension in the barostatted tilt factor must be periodic. E: Cannot use fix nvt/npt/nph with both yz dynamics and yz scaling Self-explanatory. E: Cannot use fix nvt/npt/nph with both xz dynamics and xz scaling Self-explanatory. E: Cannot use fix nvt/npt/nph with both xy dynamics and xy scaling Self-explanatory. E: Can not specify Pxy/Pxz/Pyz in fix nvt/npt/nph with non-triclinic box Only triclinic boxes can be used with off-diagonal pressure components. See the region prism command for details. E: Invalid fix nvt/npt/nph pressure settings Settings for coupled dimensions must be the same. E: Fix nvt/npt/nph damping parameters must be > 0.0 Self-explanatory. E: Cannot use fix npt and fix deform on same component of stress tensor This would be changing the same box dimension twice. E: Temperature ID for fix nvt/npt does not exist Self-explanatory. E: Pressure ID for fix npt/nph does not exist Self-explanatory. E: Fix npt/nph has tilted box too far in one step - periodic cell is too far from equilibrium state Self-explanatory. The change in the box tilt is too extreme on a short timescale. E: Could not find fix_modify temperature ID The compute ID for computing temperature does not exist. E: Fix_modify temperature ID does not compute temperature The compute ID assigned to the fix must compute temperature. W: Temperature for fix modify is not for group all The temperature compute is being used with a pressure calculation which does operate on group all, so this may be inconsistent. E: Pressure ID for fix modify does not exist Self-explanatory. E: Could not find fix_modify pressure ID The compute ID for computing pressure does not exist. E: Fix_modify pressure ID does not compute pressure The compute ID assigned to the fix must compute pressure. */ diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index e748d7f16..5f14fa23c 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -1,226 +1,228 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "stdlib.h" #include "string.h" #include "fix_wall_reflect.h" #include "atom.h" #include "comm.h" #include "update.h" #include "modify.h" #include "domain.h" #include "lattice.h" #include "input.h" #include "variable.h" #include "error.h" #include "force.h" using namespace LAMMPS_NS; using namespace FixConst; enum{XLO=0,XHI=1,YLO=2,YHI=3,ZLO=4,ZHI=5}; enum{NONE=0,EDGE,CONSTANT,VARIABLE}; /* ---------------------------------------------------------------------- */ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg < 4) error->all(FLERR,"Illegal fix wall/reflect command"); // parse args nwall = 0; int scaleflag = 1; int iarg = 3; while (iarg < narg) { if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) || (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) || (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall/reflect command"); int newwall; if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO; else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI; else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO; else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI; else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO; else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI; for (int m = 0; (m < nwall) && (m < 6); m++) if (newwall == wallwhich[m]) error->all(FLERR,"Wall defined twice in fix wall/reflect command"); wallwhich[nwall] = newwall; if (strcmp(arg[iarg+1],"EDGE") == 0) { wallstyle[nwall] = EDGE; int dim = wallwhich[nwall] / 2; int side = wallwhich[nwall] % 2; if (side == 0) coord0[nwall] = domain->boxlo[dim]; else coord0[nwall] = domain->boxhi[dim]; } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { wallstyle[nwall] = VARIABLE; int n = strlen(&arg[iarg+1][2]) + 1; varstr[nwall] = new char[n]; strcpy(varstr[nwall],&arg[iarg+1][2]); } else { wallstyle[nwall] = CONSTANT; coord0[nwall] = force->numeric(FLERR,arg[iarg+1]); } nwall++; iarg += 2; } else if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect command"); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; else error->all(FLERR,"Illegal fix wall/reflect command"); iarg += 2; } else error->all(FLERR,"Illegal fix wall/reflect command"); } // error check if (nwall == 0) error->all(FLERR,"Illegal fix wall command"); for (int m = 0; m < nwall; m++) { if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic) error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic) error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic) error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); } for (int m = 0; m < nwall; m++) if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2) error->all(FLERR, "Cannot use fix wall/reflect zlo/zhi for a 2d simulation"); // scale factors for CONSTANT and VARIABLE walls int flag = 0; for (int m = 0; m < nwall; m++) if (wallstyle[m] != EDGE) flag = 1; if (flag) { if (scaleflag) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; zscale = domain->lattice->zlattice; } else xscale = yscale = zscale = 1.0; for (int m = 0; m < nwall; m++) { if (wallstyle[m] != CONSTANT) continue; if (wallwhich[m] < YLO) coord0[m] *= xscale; else if (wallwhich[m] < ZLO) coord0[m] *= yscale; else coord0[m] *= zscale; } } // set varflag if any wall positions are variable varflag = 0; for (int m = 0; m < nwall; m++) if (wallstyle[m] == VARIABLE) varflag = 1; } /* ---------------------------------------------------------------------- */ FixWallReflect::~FixWallReflect() { + if (copymode) return; + for (int m = 0; m < nwall; m++) if (wallstyle[m] == VARIABLE) delete [] varstr[m]; } /* ---------------------------------------------------------------------- */ int FixWallReflect::setmask() { int mask = 0; mask |= POST_INTEGRATE; mask |= POST_INTEGRATE_RESPA; return mask; } /* ---------------------------------------------------------------------- */ void FixWallReflect::init() { for (int m = 0; m < nwall; m++) { if (wallstyle[m] != VARIABLE) continue; varindex[m] = input->variable->find(varstr[m]); if (varindex[m] < 0) error->all(FLERR,"Variable name for fix wall/reflect does not exist"); if (!input->variable->equalstyle(varindex[m])) error->all(FLERR,"Variable for fix wall/reflect is invalid style"); } int nrigid = 0; for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) nrigid++; if (nrigid && comm->me == 0) error->warning(FLERR,"Should not allow rigid bodies to bounce off " "relecting walls"); } /* ---------------------------------------------------------------------- */ void FixWallReflect::post_integrate() { int i,dim,side; double coord; // coord = current position of wall // evaluate variable if necessary, wrap with clear/add double **x = atom->x; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (varflag) modify->clearstep_compute(); for (int m = 0; m < nwall; m++) { if (wallstyle[m] == VARIABLE) { coord = input->variable->compute_equal(varindex[m]); if (wallwhich[m] < YLO) coord *= xscale; else if (wallwhich[m] < ZLO) coord *= yscale; else coord *= zscale; } else coord = coord0[m]; dim = wallwhich[m] / 2; side = wallwhich[m] % 2; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (side == 0) { if (x[i][dim] < coord) { x[i][dim] = coord + (coord - x[i][dim]); v[i][dim] = -v[i][dim]; } } else { if (x[i][dim] > coord) { x[i][dim] = coord - (x[i][dim] - coord); v[i][dim] = -v[i][dim]; } } } } if (varflag) modify->addstep_compute(update->ntimestep + 1); }