computing energy of a benzene crystal results in a RuntimeError

[leifjacobson]

Trying to debug some PBC issues I'm seeing with various implementations of ANI under PBC. Using an in house implementation or Torchani the energy drifts. I wanted to test with this code base but computing the energy of a benzene crystal resulted in a RuntimeError:

Traceback (most recent call last): [0/1949] File "../ani_utilities/optimize_crystal.py", line 168, in
main()
File "../ani_utilities/optimize_crystal.py", line 153, in main
energy = single_point_energy(atoms, calculator)
File "/home/jacobson/ani_utilities/utils.py", line 27, in single_point_energy
energy = molecule.get_potential_energy()/EV
File "/home/jacobson/ase_ani_venv36/lib/python3.6/site-packages/ase/atoms.py", line 664, in get_potential_e nergy
energy = self._calc.get_potential_energy(self)
File "/home/jacobson/ase_ani_venv36/lib/python3.6/site-packages/ase/calculators/calculator.py", line 505, i n get_potential_energy
energy = self.get_property('energy', atoms)
File "/home/jacobson/ase_ani_venv36/lib/python3.6/site-packages/ase/calculators/calculator.py", line 552, i n get_property
self.calculate(atoms, [name], system_changes)
File "/home/jacobson/ASE_ANI/lib/ase_interface.py", line 710, in calculate
energy, force, stddev, Fstddev = self.nc.compute_mean_props()
File "/home/jacobson/ASE_ANI/lib/ase_interface.py", line 482, in compute_mean_props
self.F[i] = nc.force().copy() RuntimeError: unidentifiable C++ exception

ERROR: CUDA throw detected! Attempting to shut down nicely! CUDA Error -- "an illegal memory access was encountered" 77 in location -- /nh/nest/u/jsmith/scratch/Gits/NeuroChem/src-atomicnnplib/cunetwork/cuannlayer_t.cu:634 in function -- m_clearDataOnDevice()

terminate called after throwing an instance of 'std::__cxx11::basic_string<char, std::char_traits, std: :allocator >' Aborted (core dumped)

The geometry I'm using is from the CCD, here is some data if you want to locate it:

data_BENZEN11 [38/1866] _audit_creation_date 2006-03-08 _database_code_depnum_ccdc_archive 'CCDC 298305' _journal_coeditor_code 'IUCr AV5045' _chemical_formula_moiety 'C6 H6' _chemical_name_systematic Benzene _journal_coden_Cambridge 622 _journal_volume 62 _journal_year 2006 _journal_page_first 94 _journal_name_full 'Acta Crystallogr.,Sect.B:Struct.Sci.'

Regards

[isayev]

Please note that ANI was not trained with PBC, so do not expect any quantitative accuracy at this point. AFAIK, benzene crystal is not orthorhombic (angles =/=90deg ). I think we don't have this functionality in the public repo. Only angles = 90 are supported now. Feel free to contact us if you are interested in this functionality.

[leifjacobson]

I understand ANI has not been trained to energies from periodic DFT and these are really just tests of the implementation not accuracy. This particular cell is orthorhombic:

_refine_ls_R_factor_gt 0.053 _refine_ls_wR_factor_gt 0.053 _diffrn_radiation_probe x-ray _symmetry_cell_setting orthorhombic _symmetry_space_group_name_H-M 'P b c a' _symmetry_Int_Tablesnumber 61 loop _symmetry_equiv_pos_site_id _symmetry_equiv_pos_as_xyz 1 x,y,z 2 1/2-x,-y,1/2+z 3 -x,1/2+y,1/2-z 4 1/2+x,1/2-y,-z 5 -x,-y,-z 6 1/2+x,y,1/2-z 7 x,1/2-y,1/2+z 8 1/2-x,1/2+y,z _cell_length_a 7.287(6) _cell_length_b 9.20(2) _cell_length_c 6.688(9) _cell_angle_alpha 90 _cell_angle_beta 90 _cell_angle_gamma 90 _cell_volume 448.366 _exptl_crystal_colour colorless _cell_formula_unitsZ 4 loop _atom_site_label _atom_site_type_symbol _atom_site_fract_x _atom_site_fract_y _atom_site_fract_z C1 C -0.0537(8) 0.1425(9) 0.0097(12) H1 H -0.085(6) 0.246(7) 0.034(8) C2 C 0.0840(7) 0.0924(10) 0.1373(10) H2 H 0.140(6) 0.156(6) 0.219(8) C3 C -0.1343(7) 0.0521(9) -0.1235(12) H3 H -0.220(6) 0.080(6) -0.204(9) C1D C 0.0537 -0.1425 -0.0097 H1D H 0.085 -0.246 -0.034 C2D C -0.084 -0.0924 -0.1373 H2D H -0.14 -0.156 -0.219 C3D C 0.1343 -0.0521 0.1235 H3D H 0.22 -0.08 0.204

[leifjacobson]

If you have an input that can be read by ASE that you have an energy and force to reproduce I'd be happy to test on my machine. Thanks for any help you can provide.

[Jussmith01]

Hi,

1) the current code only support rectangular cells. 2) if the edge length of your cell is < 2x the cutoff the code will crash with a memory error (same as your error) 3) as Oles mentioned, these potentials were not trained to work on organic crystals, so I wouldn’t expect this to work well. 4) can you define ‘the energy drifts?”. As in the energy is not conserving in NVE? If so, what precision are you running torchani at? You might want to see if switching to double helps.

[leifjacobson]

Thank you for the quick responses. I sounds like issue 2 that Justin points out is a possibility as the sidelengths of this cell are ~ 7 A whereas I believe the cutoff is ~5.5 A. I'll test with a supercell. 4. Sure, the total energy in an NVE dynamics with 0.5 fs time step is not conserved (or even close really) over the coarse of a 5 ps trajectory. Torchani is being run in double precision and the energy agrees to fairly high precision with our in house implementation which is also using double precision.

[Jussmith01]

Interesting, I’ve thoroughly tested conservation in neurochem (ASE_ANI) for gas phase systems. I also tested conservation heavily for periodic systems in recent work (https://arxiv.org/abs/1811.01914). There is some level of drift but what we saw made sense given time scale, time step, and precision of the model itself. How much drift are you seeing? Can you quantify?

[leifjacobson]

I tested a larger cubic cell of water with side length 12 A. ANI_ASE now runs without error. Tentatively the energy conservation looks better with this implementation. The disagreement in energy between this implementation and torchani is about 1.0e-4 a.u. for a single point energy evaluation on this system so something seems inconsistent. Feel free to close this issue.