Coverage for ufedmm/tests/test_respa2.py: 100%

1import numpy as np 

2import openmm 

3import pytest 

4 

5import ufedmm 

6from ufedmm.ufedmm import _standardized 

7 

8 

9def S(z): 

10 if z < 0: 

11 return 1 

12 elif z < 1: 

13 return 1 - 10 * z**3 + 15 * z**4 - 6 * z**5 

14 else: 

15 return 0 

16 

17 

18def perform_test(sigma, epsilon, charge0, charge1, rs, rc): 

19 nonbonded = openmm.NonbondedForce() 

20 nonbonded.addParticle(charge0, sigma, epsilon) 

21 nonbonded.addParticle(charge1, sigma, epsilon) 

22 nonbonded.setNonbondedMethod(nonbonded.CutoffNonPeriodic) 

23 platform = openmm.Platform.getPlatformByName("Reference") 

24 system = openmm.System() 

25 system.addParticle(1) 

26 system.addParticle(1) 

27 system.addForce(nonbonded) 

28 ufedmm.add_inner_nonbonded_force(system, rs, rc, 1) 

29 context = openmm.Context(system, openmm.CustomIntegrator(0), platform) 

30 ONE_4PI_EPS0 = 138.93545764438198 

31 for r in np.linspace(sigma, rc, 101): 

32 context.setPositions([[0, 0, 0], [r, 0, 0]]) 

33 state = context.getState(getForces=True, groups={1}) 

34 force = state.getForces()[1].x 

35 F = ( 

36 24 * epsilon * (2 * (sigma / r) ** 12 - (sigma / r) ** 6) / r 

37 + ONE_4PI_EPS0 * charge0 * charge1 / r**2 

38 ) 

39 assert force == pytest.approx(F * S((r - rs) / (rc - rs))) 

40 

41 

42def test_inner_lennard_jones(): 

43 perform_test(1.0, 1.0, 0, 0, 2.0, 2.5) 

44 

45 

46def test_inner_coulomb(): 

47 perform_test(1.0, 0.0, -1.0, 1.0, 2.0, 2.5) 

48 

49 

50def test_inner_exceptions(): 

51 model = ufedmm.AlanineDipeptideModel() 

52 nbforce = next( 

53 filter(lambda f: isinstance(f, openmm.NonbondedForce), model.system.getForces()) 

54 ) 

55 rs = 0.2 

56 rc = 0.4 

57 ufedmm.add_inner_nonbonded_force(model.system, rs, rc, 1) 

58 model.system.getForce(model.system.getNumForces() - 1).setForceGroup(3) 

59 platform = openmm.Platform.getPlatformByName("Reference") 

60 context = openmm.Context(model.system, openmm.CustomIntegrator(0), platform) 

61 context.setPositions(model.positions) 

62 forces1 = _standardized(context.getState(getForces=True, groups={3}).getForces()) 

63 forces2 = [0 * f for f in forces1] 

64 ONE_4PI_EPS0 = 138.93545764438198 

65 for index in range(nbforce.getNumExceptions()): 

66 i, j, chargeprod, sigma, epsilon = map( 

67 _standardized, nbforce.getExceptionParameters(index) 

68 ) 

69 rij = _standardized(model.positions[i] - model.positions[j]) 

70 r = np.linalg.norm(rij) 

71 z = (r - rs) / (rc - rs) 

72 F = ( 

73 S(z) 

74 * ( 

75 24 * epsilon * (2 * (sigma / r) ** 12 - (sigma / r) ** 6) / r 

76 + ONE_4PI_EPS0 * chargeprod / r**2 

77 ) 

78 * rij 

79 / r 

80 ) 

81 forces2[i] += F 

82 forces2[j] -= F 

83 for f1, f2 in zip(forces1, forces2): 

84 for i in range(3): 

85 assert f1[i] == pytest.approx(f2[i])