Issue with charges causing nonphysical geometries

[wutobias]

Description of the issue

When using LigParGen with the ACE-ALA-NME molecule, I am seeing large charges on the amide N-atoms and the H-atoms attached to these N-atoms. The charges on the N-atom are -0.997258/-1.077063 and for the H-atoms they are 0.504505/0.506935. These seem rather large in magnitude. When running a short gas-phase MD and minimizing the molecule, I am getting a structure that is ~ 8 kcal/mol lower in energy then the actual crystal structure conformation. The minimized structure looks nonphysical to me as the two amide H-atoms are pointing towards each other (see structure below). It seems to me that this issue is due to the large charges on the N/H atoms in the amide groups. My question is whether I am doing something wrong (see below for reproducing this) or if there might be bug somewhere. Any thoughts and insights on this are much appreciated.

As a site note, I am seeing the same results with CM1A-LBCC option.

Code/data to reproduce above issue

Generating force field with LigParGen:

ligpargen -i input_monomer0.pdb -cgen CM1A -c 0

Running short gas-phase MD and then minimize using openmm (v7.7.0):

def OPLS_LJ(system):
    from openmm import CustomNonbondedForce
    import numpy as np

    forces = {system.getForce(index).__class__.__name__: system.getForce(
        index) for index in range(system.getNumForces())}
    nonbonded_force = forces['NonbondedForce']
    lorentz = CustomNonbondedForce(
        '4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=sqrt(sigma1*sigma2); epsilon=sqrt(epsilon1*epsilon2)')
    lorentz.setNonbondedMethod(nonbonded_force.getNonbondedMethod())
    lorentz.addPerParticleParameter('sigma')
    lorentz.addPerParticleParameter('epsilon')
    lorentz.setCutoffDistance(nonbonded_force.getCutoffDistance())
    system.addForce(lorentz)
    LJset = {}
    for index in range(nonbonded_force.getNumParticles()):
        charge, sigma, epsilon = nonbonded_force.getParticleParameters(index)
        LJset[index] = (sigma, epsilon)
        lorentz.addParticle([sigma, epsilon])
        nonbonded_force.setParticleParameters(
            index, charge, sigma, epsilon * 0)
    for i in range(nonbonded_force.getNumExceptions()):
        (p1, p2, q, sig, eps) = nonbonded_force.getExceptionParameters(i)
        # ALL THE 12,13 and 14 interactions are EXCLUDED FROM CUSTOM NONBONDED
        # FORCE
        lorentz.addExclusion(p1, p2)
        if eps._value != 0.0:
            #print p1,p2,sig,eps
            sig14 = np.sqrt(LJset[p1][0] * LJset[p2][0])
            eps14 = np.sqrt(LJset[p1][1] * LJset[p2][1])
            nonbonded_force.setExceptionParameters(i, p1, p2, q, sig14, eps)
    return system

## Created by Leela S. Dodda for LigParGen Tutorials
## Date Aug 5, 2016
from openmm import app, KcalPerKJ
import openmm as mm
from openmm import unit as u
from sys import stdout, exit

def Minimize(simulation,iters=0):
    simulation.minimizeEnergy(maxIterations=iters)
    position = simulation.context.getState(getPositions=True).getPositions()
    energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
    app.PDBFile.writeFile(simulation.topology, position,
                          open('gasmin.pdb', 'w'))
    print ('Energy at Minima is %3.3f kcal/mol' % (energy._value * KcalPerKJ))
    return simulation

temperature = 298.15 * u.kelvin
pdb = app.PDBFile('input_monomer0.pdb')

modeller = app.Modeller(pdb.topology, pdb.positions)
forcefield = app.ForceField('input_monomer0.openmm.xml')

system = forcefield.createSystem(
    modeller.topology, nonbondedMethod=app.NoCutoff,  constraints=None)
system = OPLS_LJ(system)
integrator = mm.LangevinIntegrator(
    temperature, 1 / u.picosecond,  0.0005 * u.picoseconds)
simulation = app.Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)
simulation = Minimize(simulation,1000)

simulation.reporters.append(app.PDBReporter('gas_output.pdb', 1000))
simulation.reporters.append(app.StateDataReporter('data.txt', 1000, progress=True, temperature=True, potentialEnergy=True, density=True,totalSteps=10000,speed=True))
simulation.step(100000)

simulation = Minimize(simulation,1000)

openmm xml and pdb files: input_monomer0.pdb.gz input_monomer0.openmm.xml.gz