Implement Gaussian screening of electrostatics
It would be great to have the possibility of using screened interactions. A number of groups are working on this, e.g. http://pubs.acs.org/doi/abs/10.1021/ct5009069.
The main issue is that one extra parameter is needed per atom or atom type, which has to come from the topology in some way.
The Coulomb interaction would look like
V(r_ij) = q_i q_j erf(z_ij r_ij) / r_ij
#1 Updated by David van der Spoel almost 3 years ago
In the present form of the topology file each molecule could specify a new section with screening widths z_i. This in order to remain compatible to older topologies. E.g.:
[ atom_properties ]
; i type z_i^A z_i^B
1 1 10 12
2 1 8 6
In this way the format is ready for free energy calculations with changing z_i as well.
Since these lists would be generated normally this is not a big problem.
In order to get this into the tpr file we need to add this to the t_atoms structure, but for mdrun we need it in the mdatoms structure as well as the nbnxn_atomdata_t. That would mean minimal change to interfaces. However, a new tpr format would be needed.
#4 Updated by David van der Spoel over 2 years ago
In order to make some progress we can start by patching topology reading and introduce a new section in top files that could look like this:
[ distributed_charges ]
; Atom Type Zeta ...
1 1 5.1
2 1 6.2
where type 1 would correspond to Gaussian with one zeta.
Other types could become more complex.
#6 Updated by Mohammad Ghahremanpour about 2 years ago
David van der Spoel wrote:
The next step would be to update the SIMD kernels alternatively implement this in tables.
To use Gaussian charges, we only need the screening length, but for the Slater charges we also need to have the row number of the elements (atom) in the periodic table.
This can be determined based on the atomic number when filling the md_atoms. But, what if we want to change a Nitrogen to a Phosphor, for example, when moving the system from State A to State B. In this case, similar to mA and mB for mass, we need to have something like rowA and rowB. right?