All atom force field file generation

[saabirpetker]

Hello,

Thank you for your great software and helpful issue responses. I am having trouble generating a valid .ff file for an all-atom model of nitrile butadiene, acrylonitrile-butadiene copolymer. I would like to use parameters from https://doi.org/10.1021/acs.jpcb.6b09690 to build my copolymer systems. These parameters are in an OPLS style. Please could you assist?

Many thanks.

[fgrunewald]

Hi @saabirpetker,

It should be easy enough to generate the ff-input files so you can use polyply. I've three questions before:

1. Have you already tried generating an ff-file and could you share your progress?
2. Are these the same as default OPLS parameters or specially adopted? I didn't have time to read the entire paper
3. Are you interested in blocky copolymers or any kind of random mixing of the two units?

[saabirpetker]

Hi @saabirpetker,

Hi @fgrunewald I've only just seen this now despite living on this repo for 2 weeks - not sure how I missed it!

I've been up to solving my above problem, at least rudimentarily, to understand how I can use polyply. The tutorials were great, thanks. Your help would be much appreciated!

It should be easy enough to generate the ff-input files so you can use polyply. I've three questions before:

1. Have you already tried generating an ff-file and could you share your progress?

I'll share what I've got... It's quite basic and pretty much all manually done.
Current problems:
I didn't notice the support for dihedrals/pairs across short residues so wrote my own script to input those missing dihedral/pair parameters into my .itp files after generating them with polyply gen_params.
Also, I don't have a chemistry background, so I am sure there's a more efficient way of saying this... I am also not sure how to deal with asymmetry in residues and incorporating that linking when building a chain. For example in my case, acrylonitrile, the nitrile group carbon, "CT2" linking to the last carbon in the currently-being-built-chain, and including the associated new '; connections' interactions. I doubt it but is making two residue molecules for "forward"/"backward" facing acrylonitrile what is required?

2. Are these the same as default OPLS parameters or specially adopted? I didn't have time to read the entire paper

These are not the same, but they are not far off and I will be looking into using OPLS parameters in the future anyway, so going through that workflow would be helpful nonetheless, if that would be easier for you.

3. Are you interested in blocky copolymers or any kind of random mixing of the two units?

I'm just looking at random copolymers.

Thanks for your help,
Saabir

compressed_forcefield.zip

[fgrunewald]

@saabirpetker no worries about the delayed response. You're lucky because from tomorrow onwards I'll take a longer vacation.

I've had a look at your input files and they look pretty good. So I think you are simply missing two ingredients to deal with the dihedrals / pairs spanning more than two residues and the head-to-tail vs head-head linkage issue. Here we go:

1. I assume the linkage in the backbone is [CT1-CT2]-[CT1-CT2]-[CT1 -CT2]. Within the braces are the residues (i.e. let's consider three). To add a pair or dihedral between the CT2 of the first and the CT1 of the last residue you'd write:

[ link ]
resname "BUT|ACN"
[ atoms ]
CT2 { }
+CT1 { }
+CT2 { }
++CT1 { }
[ dihedrals ]
CT2 +CT1 +CT2 ++CT1 params {"comment": "test"}
[ pairs ]
CT2 ++CT1 1 {"comment": "test"}
[ edges ]
CT2 +CT1
+CT1 +CT2
+CT2 ++CT1

2. In order to have head-to-tail or head-to-head linkages, you could define links with edge attributes like so:

[ link ]
resname "BUT|ACN"
[ bonds ]
CT2 +CT1  1   0.1529 224262.400 {"comment": "head-to-tail"}
[ edges ]
CT2 +CT1 {"linktype": "HT"}
[ link ]
resname "BUT|ACN"
[ bonds ]
CT1 +CT1  1   0.1529 224262.400 {"comment": "head-to-head"}
[ edges ]
CT1 +CT1 {"linktype": "HH"}
[ link ]
resname "BUT|ACN"
[ bonds ]
CT2 +CT2  1   0.1529 224262.400 {"comment": "tail-to-tail"}
[ edges ]
CT2 +CT2 {"linktype": "TT"}
[ link ]
resname "BUT|ACN"
[ bonds ]
CT1 +CT2  1   0.1529 224262.400 {"comment": "tail-to-head"}
[ edges ]
CT1 +CT2 {"linktype": "TH"}

Next, you need to use gen_params with a .json sequence graph where the edges are annotated with HH or HT. This tutorial talks a bit about how to do that. To generate I wrote a little script that assigns the different linkages. I used the networkx library. The most important functions are the json export and how to set edge attributes.

You need to do this a bit manually because certain connections are not allowed. For example, you cannot have two heads connecting to the same tail at the same residue. However, you need to double-check if that this script actually does what you need.

import sys
import json
import numpy as np
import networkx as nx
from networkx.readwrite import json_graph

nmon = int(sys.argv[1])

g = nx.path_graph(nmon)
for idx, node in enumerate(g.nodes):
    resname = np.random.choice(['BUT', 'ACN'])
    g.nodes[node]['resid'] = idx+1
    g.nodes[node]['resname'] = resname

prev_connect = "none"
allowed_conect = {"H": ["TT", "TH"],
                  "T": ["HH", "HT"],
                  "none": ['HH', 'TT', 'HT', 'TH']}
for edge in nx.dfs_edges(g, source=0):
    conect = np.random.choice(allowed_conect[prev_connect])
    prev_connect = conect[1]
    g.edges[edge]['linktype'] = conect

g_json = json_graph.node_link_data(g)
with open("test.json", "w") as file_handle:
    json.dump(g_json, file_handle, indent=2)

You can call the code like so:

python gen_graph.py <number of monomers>

It will generate a residue graph with random choices of BTN and ACN and random linkage, but take into account that a head-head cannot follow another head-head. Essentially the idea is to end up with a .json input file that looks like so:

{
  "directed": false,
  "multigraph": false,
  "graph": {},
  "nodes": [
    {
      "resname": "BUT",
      "seqid": 0,
      "id": 0
    },
    {
      "resname": "ACN",
      "seqid": 0,
      "id": 1
    },
    {
      "resname": "ACN",
      "seqid": 0,
      "id": 2
    },
    {
      "resname": "ACN",
      "seqid": 0,
      "id": 3
    },
    {
      "resname": "ACN",
      "seqid": 0,
      "id": 4
    }
  ],
"links": [
    {
      "source": 0,
      "target": 1,
      "linktype": "HH"
    },
    {
      "source": 1,
      "target": 2,
      "linktype": "TH"
    },
    {
      "source": 2,
      "target": 3,
      "linktype": "TT"
    },
    {
      "source": 3,
      "target": 4,
      "linktype": "HT"

    }
  ]
}

You can now generate the polymer using the following command:

polyply gen_params -f link.ff ACN.ff BUT.ff -seqf test.json -o test.itp -name test

Enjoy! Please note that I will be back only in March and have limited access to the internet. In the meantime, @ricalessandri might help you out as well.

[saabirpetker]

Thanks @fgrunewald,

I'll give it a go. Have a good holiday!

[fgrunewald]

@saabirpetker did you have any success?

[saabirpetker]

@saabirpetker did you have any success?

@fgrunewald Thank you, yes I have. I could not quite get the ++ notation to work for me - so just wrote a Python script to do the trick, but it is a bit hacky. I should hopefully get back to you on that after my department review this July, when I may be thinking about building new structures...

I was away for quite some time, so I am still validating the NBR structures I have generated - essentially looking for changes in glass transition with changes in ACN residue wt%. This seems to be going well so far!

[saabirpetker]

Hi, I was hoping to get some help as to why polyply cannot generate structures to target densities, or even densities greater than of order 1 kgm-3, for the attached forcefields. I'm not sure what I am doing wrong here. An example .top, .json, .itp is also attached.

forcefields.zip

[fgrunewald]

Something seems wrong with the itp or force field files. The boxsize seems too small for the density requested. However, I don't have a lot of time to look into it right now. @ricalessandri perhaps you can have a look.

[saabirpetker]

forcefields+outputs.zip

Just to follow up, I tried a few things differently, but still get the same volume issues. There is attached a gro of the very low density structure that can be generated too.

[ricalessandri]

@saabirpetker let me take a look.

[ricalessandri]

@saabirpetker The problem is that you're giving the density in g/cm3 (-dens 1):

/home/saabir/Experiments/polyply-env/bin/polyply gen_coords -p ../test_auto_poly_nbr/NBR_4x250x0.4.top -o ../test_auto_poly_nbr/NBR_4x250x0.4_1kgm-3.gro -name NBR_4x250x0.4 -dens 1

while polyply needs it in kg/m3, that is, you need -dens 1000.

Take a look at:

polyply gen_coords --help

which gives:

...
  -dens DENSITY         density of system (kg/m3)
...

Please try it out and let us know.

[saabirpetker]

@ricalessandri The example above is just to show that polyply can generate a low density structure. We are both on the same page that this is not a valid structure.

When I try to generate at any higher a density, the structure is not generated, with the output remaining on "INFO - step - generating system coordinates". Even for 10 kgm-3 this step remains at 75%. Any higher and it's 0% for the 10 minutes I leave this running.

[saabirpetker]

@ricalessandri Thanks for taking a look!

[saabirpetker]

Hi, just following up to see if there are any more suggestions based on polyply getting stuck on >1kgm-3 target densities for the previous forcefield files uploaded

[fgrunewald]

Yes, @saabirpetker. Good news there is a mistake in your topology file. I think you switched the sigma epsilon columns. This mistake leads polyply to believe your atoms are 3nm in size and therefore the size of the particle in the random-walk is estimated to be about 3nm, which is non-sense of course.

Below is the more correct file. I did not check if the epsilon or sigma parameters make sense though.

[ atomtypes ]
; type name bond_type mass charge ptype sigma epsilon
BUT_CTB1  CTB1  1  12.011  -0.12  A       0.276144    3.5
BUT_CTM1  CMB1  1  12.011  -0.115  A      0.276144    3.55
BUT_CTM2  CMB2  1  12.011  -0.115  A      0.276144    3.55
BUT_CTB2  CTB2  1  12.011  -0.12  A       0.276144    3.5
BUT_HCT1  HCT1  1  1.008  0.06  A         0.12552     2.5
BUT_HCT2  HCT2  1  1.008  0.06  A         0.12552     2.5
BUT_HCM1  HCM1  1  1.008  0.115  A        0.12552     2.42
BUT_HCM2  HCM2  1  1.008  0.115  A        0.12552     2.42
BUT_HCT3  HCT3  1  1.008  0.06  A         0.12552     2.5
BUT_HCT4  HCT4  1  1.008  0.06  A         0.12552     2.5
ACN_CTB1  CTB1  1  12.011  -0.12  A       0.276144    3.5
ACN_CTN1  CTN1  1  12.011  0.04  A        0.276144    3.3
ACN_CN1  CN1  1  12.011  0.46  A          0.276144    3.3
ACN_NZ1  NZ1  1  14.007  -0.56  A         0.71128     3.2
ACN_HCT1  HCT1  1  1.008  0.06  A         0.12552     2.5
ACN_HCT2  HCT2  1  1.008  0.06  A         0.12552     2.5
ACN_HCN1  HCN1  1  1.008  0.06  A         0.06276     2.5

[saabirpetker]

I see, thank you. The issue was I had the sigma column in angstroms and not nm. So it was the right way round but the wrong scale! Now I can generate structures extremely quickly.

I plan on expanding my work to other OPLS-based parameters in the next few weeks, so will let you know of any other issues, but if you have any quick tips please let me know. Thanks again.