cabb99 / open3spn2 Goto Github PK

The following exception occurs.
Exception: CustomHbondForce: OpenCLPlatform does not support more than four exclusions per acceptor.
It may be possible to write the force in a simpler way to not have as many exclusions.

If there are more than 100K atoms in the template DNA, it causes some error

Hi, sorry if this is simple but I'm fairly new to python... I'm attempting to run the tutorial in a Jupyter notebook but in the first few lines I'm getting the following error:

UnicodeDecodeError:` 'charmap' codec can't decode byte 0x9d in position 6251: character maps to 'undefined'

`UnicodeDecodeError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_19464/3355266105.py in
2 # DNA type can be changed to 'A' or 'B'
3 seq='ATACAAAGGTGCGAGGTTTCTATGCTCCCACG'
----> 4 dna=open3SPN2.DNA.fromSequence(seq,dna_type='B_curved')
5 # Compute the topology for the DNA structure.
6 # Since the dna was generated from the sequence using X3DNA,

~\Anaconda3\envs\py37\lib\site-packages\open3SPN2\ff3SPN2.py in fromSequence(cls, sequence, dna_type, output_pdb, temp_name, compute_topology)
568 """ Initializes a DNA object from a DNA sequence """
569 self = cls()
--> 570 self.parseConfigurationFile()
571 sequence = pandas.Series([a for a in sequence], index=[('A', i) for i in range(len(sequence))])
572 # Make a possible structure

~\Anaconda3\envs\py37\lib\site-packages\open3SPN2\ff3SPN2.py in parseConfigurationFile(self, configuration_file)
157 self.configuration_file = configuration_file
158 config = configparser.ConfigParser()
--> 159 config.read(configuration_file)
160
161 # Parse all sections of the configuration file

~\Anaconda3\envs\py37\lib\configparser.py in read(self, filenames, encoding)
694 try:
695 with open(filename, encoding=encoding) as fp:
--> 696 self._read(fp, filename)
697 except OSError:
698 continue

~\Anaconda3\envs\py37\lib\configparser.py in _read(self, fp, fpname)
1012 indent_level = 0
1013 e = None # None, or an exception
-> 1014 for lineno, line in enumerate(fp, start=1):
1015 comment_start = sys.maxsize
1016 # strip inline comments

~\Anaconda3\envs\py37\lib\encodings\cp1252.py in decode(self, input, final)
21 class IncrementalDecoder(codecs.IncrementalDecoder):
22 def decode(self, input, final=False):
---> 23 return codecs.charmap_decode(input,self.errors,decoding_table)[0]
24
25 class StreamWriter(Codec,codecs.StreamWriter):

UnicodeDecodeError: 'charmap' codec can't decode byte 0x9d in position 6251: character maps to
`

I'm not getting any errors after importing all the necessary libraries... and I believe I have X3DNA correctly installed and in my env.

Desktop (please complete the following information):

• OS: Win 10
• Browser Chrome
• Version 94.0.4606.81
• Python v 3.7

Any help is greatly appreciated.

Describe the bug
After running the Protein_DNA_example.py for 10^7 steps (timestep = 2 fs), the 1lmb protein/dna complex would fully disassociate and the protein would further fold (compared with its conformation in 1lmb.pdb, where the protein extends to fit in with the major groove). To be more specific, the recognition through two major grooves disassociate fully at even 5e5 steps (~ 1 ns) and in the same time the protein folds more compactly, after which the protein only interact with the terminal of the DNA, and at around 5e6 steps the protein and DNA lose contact.

To Reproduce
Steps to reproduce the behavior:

1. open the Protein_DNA_example.py under open3spn2/example/Protein_DNA
2. edit the line 186 from simulation.step(100000) to simulation.step(10000000)
3. run the Protein_DNA_example.py

Expected behavior
As described in the OpenAWSEM/Open3SPN2 paper, I expected OpenAWSEM and Open3SPN2 make the protein/dna complex at least fluctuate around the original pdb conformation. I know that the docking is generally a hard problem, but it seems in this case the protein/dna just diffuse away from the starting conformation, indicating the pdb conformation is probably not even a local minimum? I report this issue just to make sure there is nothing I did wrong when using openawsem/open3spn2 or this is actually expected behaviour of current OpenAWSEM/Open3SPN2 (since this is actually the example provided in the paper).

Screenshots

• the start conformation
  

• the end conformation, where the protein and dna disassociate fully
  

The whole trajectory can be found here: 1lmb_traj.tar.gz

Desktop (please complete the following information):

• OS: Centos (cluster in uchicago, midway-rcc)

Dear developer,

Describe the bug
According to Eq(1) of the 3SPN2C DNA paper (S. G. Freeman et al., 2014 JCP), a cosine term of the dihedral potential is k * (1 + cos(Φ-Φ0)). However, I find that open3spn2 implements the cosine term of the dihedral potential with k * (1 - cos(Φ-Φ0)):

energy = K_periodic*(1-cs)-K_gaussian*exp(-dt_periodic^2/2/sigma^2);

Expected behavior
If we follow expression of the original 3SPN2C paper, the following expression is correct

energy = K_periodic*(1+cs)-K_gaussian*exp(-dt_periodic^2/2/sigma^2);

Could you tell me which one is correct, k * (1 + cos(Φ-Φ0)) or k * (1 - cos(Φ-Φ0)), as a potential function of the dihedral angle of the 3SPN2C DNA.

I am finding that my default versions of python, cuda, etc. computer do not meet the right dependencies for openMM or pdbfixer.

I think it would be nice if you added a quick dependencies section to the installation guide so that users know what currently works.

So far, I found that CUDA 10.2 and python 3.8 are too recent to be supported. I switched to python 3.7 with a conda install python=3.7 and am now able to install.

It might be nice to also add versions that you also notice that it works with.

The occupancy and temFactor could report error 'could not convert string to float' when there is blank.
Suggest to change like this:
occupancy=[0.0 if line[54:60].strip()=='' else float(line[54:60])],
tempFactor=[0.0 if line[60:66].strip()=='' else float(line[60:66])],

Dear Carlos,

First of all, what a nice project!
I want to simulate plasmids (dsDNA with circular topology).
However, the bonds connecting start-to-end do not seem to be connected during the simulation.

All steps followed from the example starting from PDB (and removing the protein part) work fine.
Even the simulation seems to run with constant energy and temperature.
Only after visual inspection of the trajectory I noticed the circular DNA object breaks.

I think it has to do how the pdb file deals with standard residue bonds only.
But I haven't found a way to by-pass or introduce this exception of circular topology.

Do you have any suggestions to make the code work with circular topology as well?
Of course I am happy to help!

Starting pdb and coarse grained pdb (generated with openMM) containing a plasmid of 500 bp
coarse_grained_plasmid_500.txt
atomistic_plasmid_500.txt

Describe the bug
There is always a 'segmentation fault' thrown after the simulation of 1lmb example in my local CPU platform. Although this seems harmless as we already finished the simulation, I found that if I comment out forces.update({force_name:force}) (line 179) in Protein-DNA section, this 'segmentation fault' will be thrown before the simualtion. Weirdly I actually can't reproduce the latter behaviour in my local environment but can in cluster, both platform CPU and CUDA.

To Reproduce
Steps to reproduce the behavior:

1. Go to 'open3spn2/example/Protein_DNA'
2. run jupyter nbconvert --to python Protein_DNA_example_CPU.ipynb
3. run python Protein_DNA_example_CPU.py
4. or further, comment out line 179

       174 #Add DNA-protein interaction forces
       175 for force_name in open3SPN2.protein_dna_forces:
       176     print(force_name)
       177     force = open3SPN2.protein_dna_forces[force_name](dna,protein)
       178     s.addForce(force)
       179     #forces.update({force_name: force})

5. run python Protein_DNA_example_CPU.py

Expected behavior
For the run of original Protein_DNA_example_CPU.py, we expect it to run without any errors; for commenting out the forces dictionary updating line, we expect it has no effect at all on the simulation because the dictionary forces is just used to store the names of the potentials and their force. (The reason why I comment it out is, I try to get rid of the forces entirely, but found that this segfault would be thrown earlier, and I found that the critical component is the Protein-DNA force updates. ) Practically I'll just continue to use this forces and it works, but it's an interesting (and a bit weird) bug so I am reporting it. And like I said, I can't reproduce it in local environment (even if I directly download the buggy script) but it produces bugs in the cluster.

Desktop (please complete the following information):

• OS: Manjaro (local); Scientific Linux 7.4 (cluster)

Error message
For the original run:

chain A is a DNA chain. it will be removed
chain B is a DNA chain. it will be removed
C 87
D 92
Bond
Angle
Stacking
Dihedral
BasePair
......
pap_1 term ON
No ssweight given, assume all zero
pap2
pap_2 term ON
No ssweight given, assume all zero
-899.1439809511451
TotalEnergy -899.143982 kJ/mol
Bond 327.558667 kJ/mol
Angle 973.859067 kJ/mol
.......
pap1 -0.0 kJ/mol
pap2 -0.0 kJ/mol
#"Step","Time (ps)","Potential Energy (kJ/mole)","Temperature (K)"
100,0.20000000000000015,-4280.276008936478,170.52526668552082
200,0.4000000000000003,-4134.542026299682,203.75283918870957
300,0.6000000000000004,-4036.356982095788,225.23115897642245
400,0.8000000000000006,-3966.457310441064,237.45138976781178
500,1.0000000000000007,-3889.056514398391,247.6502311417848
600,1.2000000000000008,-3818.037579398763,251.60876382973962
700,1.400000000000001,-3796.3866712439412,267.82497316002235
800,1.6000000000000012,-3821.7808458729273,282.49289463831553
900,1.8000000000000014,-3652.7479541963394,284.57682273838384
1000,2.0000000000000013,-3530.1530461379457,282.05862074578164
Bond 63.105515 kJ/mol
Angle 102.496099 kJ/mol
Stacking -446.861355 kJ/mol
Dihedral -481.204538 kJ/mol
BasePair -267.25972 kJ/mol
CrossStacking -54.91899 kJ/mol
Exclusion 1.831451 kJ/mol
Electrostatics 23.869024 kJ/mol
ExclusionProteinDNA -27.788798 kJ/mol
ElectrostaticsProteinDNA -10.197532 kJ/mol
Connectivity 1609.873732 kJ/mol
Chain 1609.873732 kJ/mol
Chi 1609.873732 kJ/mol
Excl 1609.873732 kJ/mol
rama -1587.498034 kJ/mol
rama_pro -1587.498034 kJ/mol
contact -1266.448839 kJ/mol
frag -1004.593585 kJ/mol
beta1 -184.557473 kJ/mol
beta2 -184.557473 kJ/mol
beta3 -184.557473 kJ/mol
pap1 -0.0 kJ/mol
pap2 -0.0 kJ/mol

[1]    171515 segmentation fault (core dumped)  ./Protein_DNA_example_CPU.py

For the commented run:

chain A is a DNA chain. it will be removed
chain B is a DNA chain. it will be removed
C 87
D 92
Bond
Angle
Stacking
Dihedral
BasePair
......
pap_1 term ON
No ssweight given, assume all zero
pap2
pap_2 term ON
No ssweight given, assume all zero
[2]    31762 segmentation fault  ./Protein_DNA_example_CPU.py

For this commented run, the error is thrown in the line 210:

  2 integrator = simtk.openmm.LangevinIntegrator(temperature, 1 / simtk.openmm.unit.picosecond, 2 * simtk.openmm.unit.femtoseconds)
  1 platform = simtk.openmm.Platform.getPlatformByName(platform_name)
210 simulation = simtk.openmm.app.Simulation(top,s, integrator, platform)
  1 simulation.context.setPositions(coord)
  2 energy_unit=simtk.openmm.unit.kilojoule_per_mole

I digged in a bit, and it further calls the initialization of the context, and that's the place the segfault is thrown:

pap_1 term ON
No ssweight given, assume all zero
pap2
pap_2 term ON
No ssweight given, assume all zero
> Protein_DNA_example_CPU.py(149)<module>()
-> print("Setting up the simulation...")
(Pdb) n
Setting up the simulation...
> Protein_DNA_example_CPU.py(150)<module>()
-> integrator = simtk.openmm.LangevinIntegrator(temperature, 1 / simtk.openmm.unit.picosecond, 2 * simtk.openmm.unit.femtoseconds)
(Pdb) n
> Protein_DNA_example_CPU.py(151)<module>()
-> platform = simtk.openmm.Platform.getPlatformByName(platform_name)
(Pdb) n
> Protein_DNA_example_CPU.py(152)<module>()
-> simulation = simtk.openmm.app.Simulation(top,s, integrator, platform)
(Pdb) s
--Call--
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(59)__init__()
-> def __init__(self, topology, system, integrator, platform=None, platformProperties=None, state=None):
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(82)__init__()
-> self.topology = topology
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(84)__init__()
-> if isinstance(system, string_types):
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(88)__init__()
-> self.system = system
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(90)__init__()
-> if isinstance(integrator, string_types):
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(94)__init__()
-> self.integrator = integrator
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(96)__init__()
-> self.currentStep = 0
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(98)__init__()
-> self.reporters = []
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(99)__init__()
-> if platform is None:
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(102)__init__()
-> elif platformProperties is None:
(Pdb) n
> .conda/envs/dinucl/lib/python3.6/site-packages/simtk/openmm/app/simulation.py(103)__init__()
-> self.context = mm.Context(self.system, self.integrator, platform)
(Pdb) n
[1]    201769 segmentation fault  python -m pdb Protein_DNA_example_CPU.py

It's uncertain to me why or how this bug is created, and I would really appreciate it if you can share some insights on this!

Hi Carlos,
The energy state after around 2 hundred million steps will become NaN and pops out like below. Also if I add the PDB reporter, the program also gets an error. I remembered you mentioned how to solve that, could you reply once again?

ValueError Traceback (most recent call last)
in
65 simulation.reporters.append(simtk.openmm.app.CheckpointReporter('checkpnt_1.chk', 5000))
66 #Run simulation
---> 67 simulation.step(1000000)
68

C:\Anaconda3\envs\openmm\lib\site-packages\simtk\openmm\app\simulation.py in step(self, steps)
130 def step(self, steps):
131 """Advance the simulation by integrating a specified number of time steps."""
--> 132 self._simulate(endStep=self.currentStep+steps)
133
134 def runForClockTime(self, time, checkpointFile=None, stateFile=None, checkpointInterval=None):

C:\Anaconda3\envs\openmm\lib\site-packages\simtk\openmm\app\simulation.py in _simulate(self, endStep, endTime)
233 self._generate_reports(wrapped, True)
234 if len(unwrapped) > 0:
--> 235 self._generate_reports(unwrapped, False)
236
237 def _generate_reports(self, reports, periodic):

C:\Anaconda3\envs\openmm\lib\site-packages\simtk\openmm\app\simulation.py in _generate_reports(self, reports, periodic)
252 getEnergy=getEnergy, getParameters=True, enforcePeriodicBox=periodic)
253 for reporter, next in reports:
--> 254 reporter.report(self, state)
255
256 def saveCheckpoint(self, file):

C:\Anaconda3\envs\openmm\lib\site-packages\simtk\openmm\app\statedatareporter.py in report(self, simulation, state)
193
194 # Check for errors.
--> 195 self._checkForErrors(simulation, state)
196
197 # Query for the values

C:\Anaconda3\envs\openmm\lib\site-packages\simtk\openmm\app\statedatareporter.py in _checkForErrors(self, simulation, state)
345 energy = (state.getKineticEnergy()+state.getPotentialEnergy()).value_in_unit(unit.kilojoules_per_mole)
346 if math.isnan(energy):
--> 347 raise ValueError('Energy is NaN')
348 if math.isinf(energy):
349 raise ValueError('Energy is infinite')

ValueError: Energy is NaN