No interactions between Solute and Solvent

NO_INTER_SS.py

from simtk.openmm import app,KcalPerKJ,Vec3
import simtk.openmm as mm
from simtk import unit as u
from sys import stdout,exit

#from simtk.openmm.app import *
#from simtk.openmm import *
#from simtk.unit import *
#from sys import stdout,exit
#from mdtraj.reporters import NetCDFReporter         # <-- new import from mdtraj
#import mdtraj as md

def OPLS_LJ(system):
    forces = { system.getForce(index).__class__.__name__ : system.getForce(index) for index in range(system.getNumForces()) }
    nonbonded_force = forces['NonbondedForce']
    lorentz=mm.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)
        lorentz.addExclusion(p1,p2) ## ALL THE 12,13 and 14 interactions are EXCLUDED FROM CUSTOM NONBONDED FORCE
        if eps._value != 0.0:
            sig14=u.sqrt(LJset[p1][0]*LJset[p2][0])
            eps14=u.sqrt(LJset[p1][1]*LJset[p2][1])
            nonbonded_force.setExceptionParameters(i, p1, p2, q, sig14, eps)
    return system



pdb = app.PDBFile('UNK_E1EB4D.pdb')

modeller = app.Modeller(pdb.topology, pdb.positions)
forcefield = app.ForceField('UNK_E1EB4D.xml','tip3p.xml')
modeller.addSolvent(forcefield,boxSize=Vec3(3.0, 3.0, 3.0)*u.nanometers)

system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME, nonbondedCutoff=1.0*u.nanometer, constraints=app.HBonds)
system = OPLS_LJ(system)
nb_force  = system.getForce(3)
cus_force  = system.getForce(5)
for particle1 in range(pdb.topology.getNumAtoms()):
    for particle2 in range(pdb.topology.getNumAtoms(),system.getNumParticles()):
        nb_force.addException(particle1, particle2, chargeProd=0.00, sigma=0.00, epsilon=0.00)
        cus_force.addExclusion(particle1, particle2)
integrator = mm.LangevinIntegrator(300*u.kelvin, 1/u.picosecond,  0.001*u.picoseconds)
integrator.setConstraintTolerance(0.00001)
system.addForce(mm.MonteCarloBarostat(1*u.atmospheres, 300*u.kelvin, 25))
simulation = app.Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)
## MINIMIZING 
simulation.minimizeEnergy(maxIterations=100)
position = simulation.context.getState(getPositions=True).getPositions()
app.PDBFile.writeFile(simulation.topology, position, open('boxmin.pdb', 'w'))
print('I am here')
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True,potentialEnergy=True, temperature=True,progress=True,totalSteps=100000))
simulation.step(100000)
## Setting up simulation for 10ps
TOT_STEPS = 50000000
simulation.reporters.append(app.StateDataReporter(stdout, 10, step=True,potentialEnergy=True, temperature=True,progress=True,totalSteps=TOT_STEPS))
simulation.step(TOT_STEPS)