IndexError in Neutralize Sequence

[alanrodriguez06]

Hi,

I'm trying to obtain 3D structures of a DNA oligo using the ‘free aptamer’ mode, but I’m finding a problem. 
As a test run, I set my paths and changed the sequence in main.py (all coordinates and parameters left unchanged). I’m getting an error but I can’t find a fix:   Neutralize Sequence Traceback (most recent call last):   File "/home/afrs/E2EDNA/utils.py", line 495, in findTXYZEndSoluteEnd     if (waterO in line) and (waterH in text[i + 1]) and (waterH in text[i + 2]) and (waterO in text[i + 3]):  # look for waters, atom type 349 and 350 in amoebabio18 IndexError: list index out of range   During handling of the above exception, another exception occurred:   Traceback (most recent call last):   File "main.py", line 140, in     reactionCoordinates = e2edna.runFreeAptamer() # retrieve reaction coordinates from free aptamer trajectory   File "/home/afrs/E2EDNA/e2edna.py", line 188, in runFreeAptamer     self.neutralize('sequence_soaked.xyz', 'params_combined.key')   File "/home/afrs/E2EDNA/e2edna.py", line 462, in neutralize     max_solute_index = findTXYZEndSoluteEnd(structure, self.MM_params)   File "/home/afrs/E2EDNA/utils.py", line 500, in findTXYZEndSoluteEnd     raise ValueError("No waters in structure!") ValueError: No waters in structure!  

—————

I didn’t change any parameter files and I can see   atom        349   90    O     "AMOEBA Water O"               8    15.999    2 atom        350   91    H     "AMOEBA Water H"               1     1.008    1   in the amoebabio18.prm in /params   Any thoughts on where this could come from?

Thanks!

[InfluenceFunctional]

Hi - could you send over all the params in your main.py so that we can see what the code is trying to do? Also, if the code made any .xyz files or .out files that would be very helpful for debugging.

On a semi-related note, I'll also plug our new code https://github.com/influencefunctional/opendna It is much faster, with more features and generally is more user-friendly, and is currently under active development, though you cannot use AMOEBA force-fields.

[taoliu032]

Hi there, (1) Exactly as InfluenceFunctional mentioned above, please send over the parameters, .xyz files, .out files, etc so that we can take a look. (2) To me, the error seems to indicate that the "soak" step did not successfully add water molecules to the simulation system, then the program stopped right at the beginning of the next step "neutralize". Therefore, I would pay particular attention to the "....wet.out" file in the "outfiles" directory - check if any error was reported there. (3) Meanwhile, I would also take a look at those 3 key files in your "workDir" - minimize.key, equilibrate.key, and dynamics.key - to check if the waterbox information has been filled in. If it has been filled in, then at the top of each .key file, you should see "aaxis {some_value}" instead of "aaxis XX".

PS: the error you saw meant there was no water molecules found in your simulation system, whereas the "amoebabio18.prm" is the file that provides the force field parameters for various molecules, including water.

[alanrodriguez06]

Thanks. Looking at the files and playing around I noticed there are some issues with gfortran. I'm using MMB 2.14 and tinker8. I ran the code with gnu8.3 and gnu5.4,and the error during soaking reads:

/home/afrs/tinker/bin/bin/pdbxyz: /home/afrs/Installer.2_14.Linux64/libgfortran.so.3: version 'GFORTRAN_1.4' not found (required by /home/afrs/tinker/bin/bin/pdbxyz)
Starting Fresh Run 7
Get Secondary Structure
Folding Sequence
Soak Sequence
Traceback (most recent call last):
  File "main.py", line 140, in <module>
    reactionCoordinates = e2edna.runFreeAptamer() # retrieve reaction coordinates from free aptamer trajectory
  File "/home/afrs/E2EDNA/e2edna.py", line 178, in runFreeAptamer
    copyfile('sequence.xyz', 'sequence_pre_soak.xyz')
  File "/home/afrs/.conda/envs/MD/lib/python3.6/shutil.py", line 120, in copyfile
    with open(src, 'rb') as fsrc:
FileNotFoundError: [Errno 2] No such file or directory: 'sequence.xyz'

The files generated are attached in Folder 1. Folder-1 MD.zip

Then I tried linking to the library I already had as suggested by others (https://stackoverflow.com/questions/9628273/libgfortran-version-gfortran-1-4-not-found) :

ln -sf /usr/lib64/libgfortran.so.3 /home/afrs/Installer.2_14.Linux64/libgfortran.so.3

And got this:

At line 654 of file /home/afrs/tinker/source/xyzedit.f (unit = 5, file = 'stdin')
Fortran runtime error: End of file
At line 993 of file /home/afrs/tinker/source/xyzedit.f (unit = 5, file = 'stdin')
Fortran runtime error: End of file
Starting Fresh Run 10
Get Secondary Structure
Folding Sequence
Soak Sequence
Neutralize Sequence
Traceback (most recent call last):
  File "/home/afrs/E2EDNA/utils.py", line 495, in findTXYZEndSoluteEnd
    if (waterO in line) and (waterH in text[i + 1]) and (waterH in text[i + 2]) and (waterO in text[i + 3]):  # look for waters, atom type 349 and 350 in amoebabio18
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "main.py", line 140, in <module>
    reactionCoordinates = e2edna.runFreeAptamer() # retrieve reaction coordinates from free aptamer trajectory
  File "/home/afrs/E2EDNA/e2edna.py", line 188, in runFreeAptamer
    self.neutralize('sequence_soaked.xyz', 'params_combined.key')
  File "/home/afrs/E2EDNA/e2edna.py", line 462, in neutralize
    max_solute_index = findTXYZEndSoluteEnd(structure, self.MM_params)
  File "/home/afrs/E2EDNA/utils.py", line 500, in findTXYZEndSoluteEnd
    raise ValueError("No waters in structure!")
ValueError: No waters in structure!

The files generated are attached in Folder 2. Folder-2 MD.zip

The problem seems to be MMB can't use xyzedit to create the coordinates during soaking. Actually, the wet.out file was not generated when using the gfortran linked from my system.

[taoliu032]

Please double check if the two zip files indeed belong to the two runs, respectively. Thanks.

I found them quite perplexing: e.g., in "Folder-1 MD", there is "make_sequence.xyz_water_box.out", which should be generated when the method "soak" was called on line 179: self.soak('sequence.xyz', 'params_combined.key') However, the program had stopped at line 178 copyfile('sequence.xyz', 'sequence_pre_soak.xyz') according to the error. Then how was the "make_sequence.xyz_water_box.out" generated? Similarly, I would expect more outputs from your second example than those included in "Folder-2 MD". I am guessing that these two zip files might not come from the two examples in the current order.

[alanrodriguez06]

Sorry for the confusion. Here I attach 2 folders, each with the different errors and their corresponding run files. Thank you!

Error1.zip Error2.zip

[taoliu032]

In Error2 folder: as the error file says, the tinker program "pdbxyz" is missing, so that even if there was "sequence.pdb" after foldSequence(), the .pdb file was not converted to .xyz file. Consistently, outfiles/pdbconversion.out is empty, indicating the conversion from .pdb to .xzy never happened.

In Error1 folder: also as the error file indicates, the tinker program "xyzedit" is problematic, so the soaking step is not completed at all - look at the "sequence_soaked.xyz".

Consider contacting the administrator(s) of the cluster that you used to check if the tinker programs were properly installed. Meanwhile, please take a minute to read through our code which should help with your debugging.

PS: the "params" folder you attached is the force field file for AMOEBA18, whereas the "params" we meant earlier were referred to the parameters in the main.py, which were specified for your simulation runs.

[alanrodriguez06]

I reinstalled all the packages. I carried out some tests using the preloaded tests on tinker8, and the programs seem to be working fine. However, I still get this error:

At line 654 of file /home/afrs/tinker/source/xyzedit.f (unit = 5, file = 'stdin')
Fortran runtime error: End of file
At line 993 of file /home/afrs/tinker/source/xyzedit.f (unit = 5, file = 'stdin')
Fortran runtime error: End of file
Starting Fresh Run 49
Get Secondary Structure
Folding Sequence
Soak Sequence
Neutralize Sequence
Traceback (most recent call last):
  File "/home/afrs/E2EDNA/utils.py", line 495, in findTXYZEndSoluteEnd
    if (waterO in line) and (waterH in text[i + 1]) and (waterH in text[i + 2]) and (waterO in text[i + 3]):  # look for waters, atom type 349 and 350 in amoebabio18
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "main.py", line 140, in <module>
    reactionCoordinates = e2edna.runFreeAptamer() # retrieve reaction coordinates from free aptamer trajectory
  File "/home/afrs/E2EDNA/e2edna.py", line 188, in runFreeAptamer
    self.neutralize('sequence_soaked.xyz', 'params_combined.key')
  File "/home/afrs/E2EDNA/e2edna.py", line 462, in neutralize
    max_solute_index = findTXYZEndSoluteEnd(structure, self.MM_params)
  File "/home/afrs/E2EDNA/utils.py", line 500, in findTXYZEndSoluteEnd
    raise ValueError("No waters in structure!")
ValueError: No waters in structure!

Fortran runtime error: End of file occurs when the file cannot be found. The sequence.xyz is formed at the end of fold, but there is a problem during soaking. I deleted the 'waterbox.in' reset, and indeed the file is edited with the x,y,z ranges. I suspect it may have something to do with the parameters file. I'm not using a ligand, so I deleted the if condition in main.py:

if params['target'] == 'UTP-4':

and just set all the parameters. I deleted the info from the ligand in the minimize, dynamics and equilibirum keys. But I'm not sure what would go into the params_combined.key, as they are input in the soak function. Also, I see at the end of the soaking step these 3 keys are reset, but they are not input of the soak function. Can you clarify on this part? Thanks! Error.zip

[InfluenceFunctional]

Sorry to keep asking but if you could upload all the .key, .in and .out files for the run(s) we're talking about? That would help a lot, as contain the debug outputs. As a general question, are you able to run Tinker utilities such as xyzedit, pdbxyz, etc. from the command line?

As to the UTP parameters in the .key files - Indeed we noticed that you had them in. They shouldn't cause problems unless you have a conflicting ligand somewhere in the simulation, but you should also have no issues if you delete them.

[alanrodriguez06]

Thank you. I attach the following:

1) I suspected the problem came from the parameter (and not the tinker programs) because the file sequence.xyz is generated (which uses the program pdbxyz). As well, I carried out a few tests included in the tinker8 files. For example, in command line in /home/afrs/runs/pdbxyztest/ (which contains the pdb and key files):

chmod +x /home/afrs/tinker/bin/pdbxyz /home/afrs/tinker/bin/pdbxyz sequence -k pdbxyztest.key > pdbconversion.out

pdbxyztest.zip

2) Running main.py gets the same error of Fortran runtime error: End of file which is because it cannot open the file, likely it doesn't exist. Here, I deleted the condition if params['target'] == 'UTP-4': in main.py but still included the path to the ligand otherwise there's an error in e2edna.py when trying to copy the file (line 38). Note that when run in the foreground I have to quit the program when it is stuck in 'Soak sequence', otherwise it just wouldn't move. When run in the background, it quits automatically during soaking showing the error I included in the previous message. Error1.zip

I see the soak function takes params_combined.key and reset the minimize, dynamics and, equilibrium keys at the end. I don’t quite understand how to configure params_combined.key so the latter is possible.

[InfluenceFunctional]

For params_combined, simply deleting the UTP-4 info, as it looks like you have done, should be fine. Again, if no atoms with the relevant parameters are found in the sim, these should not play a role.

The fact that it hangs when run in the foreground, plus the error outputs you have sent over strongly suggest that this line in e2edna.py is hanging/failing, which often happens when Tinker receives arguments which it doesn't know how to parse

os.system(self.params['xyzedit path'] + ' ' + structure + ' -k ' + key + ' 20 ' + 'water.xyz_2 > outfiles/' + structure + 'wet.out') # place structure in box of waters

To debug this, I would run xyzedit interactively on the command line and see what kind of errors it spits out. Usually these would show up in the .out file, so we may be seeing a new kind of error from Tinker here (at least in our experience). The way one parses the command line arguments can sometimes be finnicky, perhaps even between different versions of Tinker.

Something like xyzedit sequence.xyz -k params_combined.key with some sequence.xyz and keyfile in the local directory (and amoebabio18.prm in the location specified in the keyfile) then work interactively in the prompt to add a solvent box, and see if it works. If there is an error that will be a clue to fixing the issue. If it works interactively, but not via the os.system() call, we may have to engineer a new way of passing arguments to this function for your version of Tinker.

[alanrodriguez06]

Thanks for pointing this out, I didn't consider that. In fact, I can use those programs in the command line with the required files in the local directory. However, I found that the flags in e2edna.py do not match exactly what's required in each option of xyzedit. For example:

        os.system(self.params['xyzedit path'] + ' ' + structure + ' -k ' + key + ' 12 0 > outfiles/' + structure + '_centered.out') # make sure structure is centered
        removetwo(structure)
        os.system(self.params['xyzedit path'] + ' ' + structure + ' -k ' + key + ' 20 ' + 'water.xyz_2 > outfiles/' + structure + 'wet.out') # place structure in box of waters
        removetwo(structure)

I used the package hanging_threads and saw the program hangs just before centering (in foreground). When I checked the xyzedit.f (attached here), I see that option (12) is Translate All Atoms by an X,Y,Z-Vector', and option 20 is Make Rhombic Dodecahedron from Cubic Box'. Maybe we have different xyzedit files? My guess is I need to use Translate Center of Mass to the Origin' (option 13 in my file) and Soak Current Molecule in Box of Solvent' (option 23). As well, the addition of ions in the neutralize function, my guess is it should be the option Place Monoatomic Ions around a Solute' (24 in my file but 21 in yours?).

xyzedit.f.zip

I made these changes and the soak step is completed. Now I get an error during neutralize:

Soak Sequence
Neutralize Sequence
Traceback (most recent call last):
  File "main.py", line 137, in <module>
    reactionCoordinates = e2edna.runFreeAptamer() # retrieve reaction coordinates from free aptamer trajectory
  File "/home/afrs/E2EDNA/e2edna.py", line 190, in runFreeAptamer
    self.neutralize('sequence_soaked.xyz', 'params_combined.key')
  File "/home/afrs/E2EDNA/e2edna.py", line 464, in neutralize
    max_solute_index = findTXYZEndSoluteEnd(structure, self.MM_params)
  File "/home/afrs/E2EDNA/utils.py", line 502, in findTXYZEndSoluteEnd
    solLine = max_solute_line.split(' ')
UnboundLocalError: local variable 'max_solute_line' referenced before assignment

I tried making the solLine a global variable in utils.py but didn't work. Any suggestions? Thank you very much! Error2.zip

[InfluenceFunctional]

Ah - it seems indeed that we are working with different versions of Tinker. This is surprising since we have worked with several versions of Tinker 8 and Tinker 9, and never had this issue, in any case it this is a good step towards debugging.

Thanks for enclosing your zipped directory. It looks here that the issue is not in neutralization but still in the soaking. If you look at sequence_soaked.zyz, there are no waters added, just a single water oxygen at coordinates [0,0,0]. It looks to me that the water box was created improperly as well (by another xyzedit call)

os.system(self.params['xyzedit path'] +' < waterbox.in > outfiles/make_' + structure + '_water_box.out') # make box of waters Perhaps waterbox.in is bad here, or as before the interactive controls are wrong.