AbramsGroup / HTPolyNet

Development of High-Throughput Polymer Network Atomistic Simulation
https://abramsgroup.github.io/HTPolyNet/
MIT License
15 stars 7 forks source link

Polymerizations in water. #4

Closed ggolde closed 5 months ago

ggolde commented 6 months ago

First off, thank you for all the effort you put in creating and maintaining this package. I am using HTPolyNet to investigate the properties of poly(acrylamide-co-N,N'-methylenebisacrylamide) hydrogel. I have successfully generated the crosslinked polymer with no water present in the system, but since I am investigating swelling behavior, I am trying to simulate the polymerization in solution.

I have tried adding water directly to the constituents directive in the .yaml config file as if it were a monomer with no reactions, but receive the error:

File "/home/griffin/anaconda3/envs/mol-env/lib/python3.12/site-packages/HTPolyNet/topology.py", line 74, in typedata return int(s) ^^^^^^ ValueError: invalid literal for int() with base 10: '#ifdef'

I believe the error is happening when HTPolyNet is reading back in the H2O.top file it generates in proj-*/molecules/parameterized/ with the Topology class method read_top(). The generated H2O.top file contains the line "#ifdef FLEXIBLE" under the "[ atoms ]" directive that read_top() tries to convert to type int because it thinks the line is a row entry in under the [ atoms ] directive. Would modifying line 193 in Topology.py work to avoid this error (shown below)? Would it cause other errors?

Change line 193 in Topology.py from: "contentlines=[l for l in s[1:] if not (len(l.strip())==0 or l.strip().startswith(';'))]" to something like: "contentlines=[l for l in s[1:] if not (len(l.strip())==0 or ( l.strip().startswith(';') or l.strip().startswith('#') ) )]"

I am still an undergraduate and inexperienced in GROMACS. Are there any other issues you foresee by trying to simulate polymerization in solution? Would it be better to use "gmx solvate" after densification and before CURE? Thank you again for all your hard work.

My log files as well as my input H2O.mol2 file and the parameterized H2O.top file are attached (as .txt). console.log diagnostics.log H2O_mol2.txt H2O_top.txt

cameronabrams commented 6 months ago

Hi Griffin,

Thanks for bringing this to my attention; I've never tried something like this. It looks to me like gromacs requires specific water parameters depending upon whether holonomic bondlength restraints are being used (the settle algorithm). If you modified the top file so that the #-lines are removed including the entire # else branch with the [settles] stanza, this would fix the problem, but then gromacs would run with flexible water (this is not bad, just inefficient). Doing the complementary thing will not work since the topology parser in HTPolyNet does not know about [settles] stanzas.

Your idea of using gmx solvate is good; it is out of the workflow of HTPolyNet, but should be ok. You could do a build through densification setting the pressure during densification to something like 10^-5, and then gmx solvate the output of that, and then do your own equilibrations after that.

I'll add this problem to the current buglist/future features list and hopefully I can find some time to get around to this soon; let me know if the gmx solvate thing works.

Cheers, Cameron

On Mon, May 20, 2024 at 4:19 PM Griffin Golde @.***> wrote:

First off, thank you for all the effort you put in creating and maintaining this package. I am using HTPolyNet to investigate the properties of poly(acrylamide-co-N,N'-methylenebisacrylamide) hydrogel. I have successfully generated the crosslinked polymer with no water present in the system, but since I am investigating swelling behavior, I am trying to simulate the polymerization in solution.

I have tried adding water directly to the constituents directive in the .yaml config file as if it were a monomer with no reactions, but receive the error:

File "/home/griffin/anaconda3/envs/mol-env/lib/python3.12/site-packages/HTPolyNet/topology.py", line 74, in typedata return int(s) ^^^^^^ ValueError: invalid literal for int() with base 10: '#ifdef'

I believe the error is happening when HTPolyNet is reading back in the H2O.top file it generates in proj-*/molecules/parameterized/ with the Topology class method read_top(). The generated H2O.top file contains the line "#ifdef FLEXIBLE" under the "[ atoms ]" directive that read_top() tries to convert to type int because it thinks the line is a row entry in under the [ atoms ] directive. Would modifying line 193 in Topology.py work to avoid this error (shown below)? Would it cause other errors?

Change line 193 in Topology.py from: "contentlines=[l for l in s[1:] if not (len(l.strip())==0 or l.strip().startswith(';'))]" to something like: "contentlines=[l for l in s[1:] if not (len(l.strip())==0 or ( l.strip().startswith(';') or l.strip().startswith('#') ) )]"

I am still an undergraduate and inexperienced in GROMACS. Are there any other issues you foresee by trying to simulate polymerization in solution? Would it be better to use "gmx solvate" after densification and before CURE? Thank you again for all your hard work.

My log files as well as my input H2O.mol2 file and the parameterized H2O.top file are attached (as .txt). console.log https://github.com/AbramsGroup/HTPolyNet/files/15381450/console.log diagnostics.log https://github.com/AbramsGroup/HTPolyNet/files/15381451/diagnostics.log H2O_mol2.txt https://github.com/AbramsGroup/HTPolyNet/files/15381452/H2O_mol2.txt H2O_top.txt https://github.com/AbramsGroup/HTPolyNet/files/15381453/H2O_top.txt

— Reply to this email directly, view it on GitHub https://github.com/AbramsGroup/HTPolyNet/issues/4, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFHYWXIXLTIY5K7CWJIQ2TZDJLD3AVCNFSM6AAAAABIAKIBK6VHI2DSMVQWIX3LMV43ASLTON2WKOZSGMYDMNZSGIYTAOA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

ggolde commented 5 months ago

Hi Cameron,

After some fiddling, I was successfully able to run a polymerization with flexible tip3p water (named SOL) with the attached parameterized inputs. This only works if the SOL constituent is below the other monomer constituents in the "system".yaml file. If the SOL constituent was placed before/above the monomers, the atom indexes for the [ pairs ] and [ dihedrals ] directives in the init.top file were not updated leading to "No default Proper Dih. types" GROMACS errors. I have not yet located the source of this bug. I tried to resolve it by putting empty [ pairs ] and [ dihedrals ] directive in SOL.top and an empty [ dihedraltypes ] directive in SOL.itp. This successfully caused the atom indexes in init.top to be updated, but the entries in [ pairs ], [ dihedrals ], and [ dihedraltypes ] printed in floats, causing GROMACS errors.

I am going to play around with a fork of HTPolyNet to try to resolve these bugs and make it recognize the [ settles ] and [ exclusion ] directives so simulations with rigid water can be performed.

I also tried the approach with "gmx solvate," but the gmx solvate command does not correctly update the topology file with new atoms unless a "tip3p.itp" file is included but then there are errors with multiple [ atomtypes ] and [ bondtypes ] directives.

Thanks, Griffin Golde

SOL_gro.txt SOL_grx.txt SOL_itp.txt SOL_mol2.txt SOL_top.txt SOL_tpx.txt