Merge pull request #856 from akohlmey/fix_srp_bugfix

Bugfix for SRP with triclinic cells
This commit is contained in:
Steve Plimpton 2018-03-29 08:10:55 -06:00 committed by GitHub
commit 2a930dac42
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 15 additions and 6569 deletions

File diff suppressed because it is too large Load Diff

View File

@ -1,32 +0,0 @@
# Solvated 5-mer peptide
# Demonstrating temper_npt
units real
atom_style full
pair_style lj/charmm/coul/long 8.0 10.0 10.0
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic
kspace_style pppm 0.0001
read_data data.peptide
neighbor 2.0 bin
neigh_modify delay 5
timestep 2.0
thermo_style custom step temp epair emol etotal press density
thermo 50
variable temper_T world 275 280 285 290 295 300 305 310
variable rep world 0 1 2 3 4 5 6 7
fix myfix all npt temp ${temper_T} ${temper_T} 100.0 iso 1 1 1000
run 500
temper_npt 2000 100 ${temper_T} myfix 0 58728 1
fix 2 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
group peptide type <= 12

View File

@ -264,12 +264,21 @@ void FixSRP::setup_pre_force(int zz)
rmax = sqrt(rsqmax); rmax = sqrt(rsqmax);
double cutneighmax_srp = neighbor->cutneighmax + 0.51*rmax; double cutneighmax_srp = neighbor->cutneighmax + 0.51*rmax;
// find smallest cutghost double length0,length1,length2;
double cutghostmin = comm->cutghost[0]; if (domain->triclinic) {
if (cutghostmin > comm->cutghost[1]) double *h_inv = domain->h_inv;
cutghostmin = comm->cutghost[1]; length0 = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]);
if (cutghostmin > comm->cutghost[2]) length1 = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]);
cutghostmin = comm->cutghost[2]; length2 = h_inv[2];
} else length0 = length1 = length2 = 1.0;
// find smallest cutghost.
// comm->cutghost is stored in fractional coordinates for triclinic
double cutghostmin = comm->cutghost[0]/length0;
if (cutghostmin > comm->cutghost[1]/length1)
cutghostmin = comm->cutghost[1]/length1;
if (cutghostmin > comm->cutghost[2]/length2)
cutghostmin = comm->cutghost[2]/length2;
// stop if cutghost is insufficient // stop if cutghost is insufficient
if (cutneighmax_srp > cutghostmin){ if (cutneighmax_srp > cutghostmin){