Skip to content

Commit

Permalink
Fixed some issues in sc forces when forces are zero
Browse files Browse the repository at this point in the history
  • Loading branch information
@ceriottm
ceriottm committed Nov 24, 2023
1 parent dd6a884 commit e31f872
Show file tree
Hide file tree
Showing 4 changed files with 38 additions and 34 deletions.
4 changes: 2 additions & 2 deletions examples/profiling/remd_scpimd_gas/init.xyz
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
8
# CELL(abcABC): 50.0 50.0 50.0 90.00000 90.00000 90.00000
H 0 0 0
# CELL(abcABC): 10.0 10.0 10.0 90.00000 90.00000 90.00000
H 0 0 1
H 0 0 2
H 0 0 3
H 0 0 4
H 0 0 5
H 0 0 6
H 0 0 7
H 0 0 8
26 changes: 13 additions & 13 deletions examples/profiling/remd_scpimd_gas/input.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,12 +17,12 @@
</prng>
<system_template>
<labels> [ PREFIX, FILENAME, TEMPERATURE, PRESSURE ] </labels>
<instance> [ REF00, init.xyz, 300, 1 ] </instance>
<instance> [ REF01, init.xyz, 310, 1 ] </instance>
<instance> [ REF02, init.xyz, 320, 1 ] </instance>
<instance> [ REF10, init.xyz, 300, 10 ] </instance>
<instance> [ REF11, init.xyz, 310, 10 ] </instance>
<instance> [ REF12, init.xyz, 320, 10 ] </instance>
<instance> [ REF00, init.xyz, 1, 1 ] </instance>
<instance> [ REF01, init.xyz, 2, 1 ] </instance>
<instance> [ REF02, init.xyz, 5, 1 ] </instance>
<instance> [ REF10, init.xyz, 1, 10 ] </instance>
<instance> [ REF11, init.xyz, 2, 10 ] </instance>
<instance> [ REF12, init.xyz, 5, 10 ] </instance>
<template>
<system prefix="PREFIX">
<initialize nbeads="4">
Expand All @@ -38,21 +38,21 @@
<motion mode="dynamics">
<dynamics mode="scnpt">
<barostat mode='sc-isotropic'>
<tau units='femtosecond'> 200</tau>
<tau> 10</tau>
<thermostat mode='langevin'>
<tau units='femtosecond'> 100</tau>
<tau> 1</tau>
</thermostat>
<h0> [ 32.84024, 0, 0, 0, 32.84024, 0, 0, 0, 32.84024 ]</h0>
<h0> [ 10, 0, 0, 0, 10, 0, 0, 0, 10 ]</h0>
</barostat>
<thermostat mode='langevin'>
<tau units='femtosecond'> 100 </tau>
<tau> 1 </tau>
</thermostat>
<timestep units="femtosecond"> 1.0 </timestep>
<timestep> 0.01 </timestep>
</dynamics>
</motion>
<ensemble>
<temperature units="kelvin"> TEMPERATURE </temperature>
<pressure units="megapascal"> PRESSURE </pressure>
<temperature> TEMPERATURE </temperature>
<pressure> PRESSURE </pressure>
<bias> <force forcefield="dummy.2"/> </bias>
<bias_weights> [0.0] </bias_weights>
</ensemble>
Expand Down
26 changes: 13 additions & 13 deletions examples/profiling/remd_scpimd_noff/input.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,12 @@
</prng>
<system_template>
<labels> [ PREFIX, FILENAME, TEMPERATURE, PRESSURE ] </labels>
<instance> [ REF00, init.xyz, 300, 1 ] </instance>
<instance> [ REF01, init.xyz, 310, 1 ] </instance>
<instance> [ REF02, init.xyz, 320, 1 ] </instance>
<instance> [ REF10, init.xyz, 300, 10 ] </instance>
<instance> [ REF11, init.xyz, 310, 10 ] </instance>
<instance> [ REF12, init.xyz, 320, 10 ] </instance>
<instance> [ REF00, init.xyz, 1, 1 ] </instance>
<instance> [ REF01, init.xyz, 2, 1 ] </instance>
<instance> [ REF02, init.xyz, 5, 1 ] </instance>
<instance> [ REF10, init.xyz, 1, 10 ] </instance>
<instance> [ REF11, init.xyz, 2, 10 ] </instance>
<instance> [ REF12, init.xyz, 5, 10 ] </instance>
<template>
<system prefix="PREFIX">
<initialize nbeads="4">
Expand All @@ -36,21 +36,21 @@
<motion mode="dynamics">
<dynamics mode="scnpt">
<barostat mode='sc-isotropic'>
<tau units='femtosecond'> 200</tau>
<tau> 10</tau>
<thermostat mode='langevin'>
<tau units='femtosecond'> 100</tau>
<tau> 1</tau>
</thermostat>
<h0> [ 32.84024, 0, 0, 0, 32.84024, 0, 0, 0, 32.84024 ]</h0>
<h0> [ 10, 0, 0, 0, 10, 0, 0, 0, 10 ]</h0>
</barostat>
<thermostat mode='langevin'>
<tau units='femtosecond'> 100 </tau>
<tau> 1 </tau>
</thermostat>
<timestep units="femtosecond"> 1.0 </timestep>
<timestep> 0.01 </timestep>
</dynamics>
</motion>
<ensemble>
<temperature units="kelvin"> TEMPERATURE </temperature>
<pressure units="megapascal"> PRESSURE </pressure>
<temperature> TEMPERATURE </temperature>
<pressure> PRESSURE </pressure>
<bias> <force forcefield="dummy.2"/> </bias>
<bias_weights> [0.0] </bias_weights>
</ensemble>
Expand Down
16 changes: 10 additions & 6 deletions ipi/engine/forces.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1120,10 +1120,14 @@ def forcesvirs_4th_order(self, index):
# calculates the finite displacement.
fbase = dstrip(self.f)
eps = self.mforces[index].epsilon
delta = np.abs(eps) / np.sqrt(
(fbase / self.beads.m3 * fbase / self.beads.m3).sum()
/ (self.nbeads * self.natoms)
)
foverm = np.sqrt((fbase / self.beads.m3 * fbase / self.beads.m3).sum() / (self.nbeads * self.natoms))
if np.abs(foverm)>1e-20:
delta = np.abs(eps) / np.sqrt(
(fbase / self.beads.m3 * fbase / self.beads.m3).sum()
/ (self.nbeads * self.natoms)
)
else: # defaults to eps if otherwise we'd get an 1/0
delta = np.abs(eps)
dq = delta * fbase / self.beads.m3

# stores the force component.
Expand All @@ -1141,11 +1145,11 @@ def forcesvirs_4th_order(self, index):
# This makes the FD incorrect. A fix could be to ALWAYS compute the odd and the even
# beads on two different ring polymers of P / 2 beads, but centered difference + RPC
# seems to be a neater solution. Anyway, the cost of a CNT-DIFF in comparison to a FWD-DIFF
# is marginal in when RPC + MTS is used.
# is marginal when RPC + MTS is used.

if self.mforces[index].nbeads != self.nbeads:
warning(
"ERROR: high order PIMD + RPC works with a centered finite difference only! (Uness you find an elegant solution :))"
"ERROR: high order PIMD + RPC works with a centered finite difference only! (Unless you find an elegant solution :))"
)
exit()

Expand Down

0 comments on commit e31f872

Please sign in to comment.