Running the Build¶
Now, in our working directory my_pms/2-polymethylstyrene/, if you haven’t already, go ahead and launch via ./run.sh.
$ ./run.sh &
As the htpolynet run command inside run.sh indicates, standard output is being redirected from the console to the file console.log.
Console output¶
In this section, let’s go section by section through the console log console.log to understand what HTPolyNet is doing. The first section just shows a banner and indicates where the library is and where the project directory is, then reports on the available Ambertools software, finally ending by stating the name of the configuration file.
INFO>
INFO> HTPolyNet 1.0.8
INFO> https://abramsgroup.github.io/HTPolyNet/
INFO>
INFO> Ming Huang
INFO> mh3429@dragons.drexel.edu
INFO>
INFO> Cameron F. Abrams
INFO> cfa22@drexel.edu
INFO>
INFO> Supported in part by Grants W911NF-17-2-0227
INFO> and W911NF-12-R-0011 from the US Army Research Lab
INFO>
INFO> Please cite the HTPolyNet paper:
INFO>
INFO> Ming Huang and Cameron F. Abrams, HTPolyNet: A general
INFO> system generator for all-atom molecular simulations of
INFO> amorphous crosslinked polymers, SoftwareX, vol. 21,
INFO> pp. 101303, 2023 (doi:10.1016/j.softx.2022.101303)
INFO>
INFO> ******************** HTPolyNet runtime begins *********************
INFO> New project in /home/cfa/htpolynet-tutorials/1.0.8/2-polymethylstyrene/proj-0
INFO> *************************** Ambertools: ***************************
INFO> ******************** antechamber (ver. 22.0) *********************
INFO> ******************** tleap (ver. 22.0) *********************
INFO> ******************** parmchk2 (ver. 22.0) *********************
INFO> Configuration: pMSTY.yaml
Next comes the monomer and oligomer template parameterizations. HTPolyNet first handles the molecules that are explicitly named in the configuration. Only EMB.mol2 exists; the other two are generated automatically using their respective reaction directives. The three molecules implied by the concept of “chaining” are also generated and parameterized.
INFO> *********** Templates in proj-0/molecules/parameterized ***********
INFO> 3 molecules detected in pMSTY.yaml
INFO> explicit: 3
INFO> implied by stereochemistry: 0
INFO> implied by symmetry: 0
INFO> AmberTools> generating GAFF parameters from EMB.mol2
INFO> EMB: 120.19 g/mol
INFO> AmberTools> generating GAFF parameters from EMB~C1-C2~EMB.mol2
INFO> EMB~C1-C2~EMB: 238.36 g/mol
INFO> AmberTools> generating GAFF parameters from EMBCC.mol2
INFO> EMBCC: 118.17 g/mol
INFO> 3 molecules implied by chaining
INFO> AmberTools> generating GAFF parameters from EMB~C1=C2~EMB~C1-C2~EMB.mol2
INFO> EMB~C1=C2~EMB~C1-C2~EMB: 356.53 g/mol
INFO> AmberTools> generating GAFF parameters from EMB~C1-C2~EMB~C1=C2~EMB.mol2
INFO> EMB~C1-C2~EMB~C1=C2~EMB: 356.53 g/mol
INFO> AmberTools> generating GAFF parameters from EMB~C1-C2~EMB~C1=C2~EMB~C1-C2~EMB.mol2
INFO> EMB~C1-C2~EMB~C1=C2~EMB~C1-C2~EMB: 474.70 g/mol
INFO> Generated 6 molecule templates
INFO> Initial composition is EMB 127
INFO> 100% conversion is 127 bonds
If we look in proj-0/molecules/parameterized we’ll see the gro, itp, top, grx, and tpx files for each molecule. The first three are Gromacs-specific. The grx file contains “extended attributes” of each atom that HTPolyNet uses internally and are not needed for Gromacs. The tpx file indentifies topological features of interest for HTPolyNet but not necessary for Gromacs; for now, this is limited to identification of ring systems via atom indices.
Next comes initialization of the system topology and coordinates. Here, using the constituents directive, HTPolyNet generates a full system topology and simulation box. The box is filled according to the initial_density subdirective of the densification directive.
INFO> ************** Initialization in proj-0/systems/init **************
INFO> Topology "init.top" in proj-0/systems/init
INFO> Initial density: 300.0 kg/m^3
INFO> Total mass: 2.535e-23 kg
INFO> Box aspect ratio: 1.0 x 1.0 x 1.0
INFO> Initial box side lengths: 4.388 nm x 4.388 nm x 4.388 nm
INFO> Coordinates "init.gro" in proj-0/systems/init
INFO> Extended attributes "init.grx" in proj-0/systems/init
Next comes a report of the densification of the system. Note that because of the stochastic nature of initializing and equilibrating, your numbers for density and box side lengths may be a little different from those shown below.
INFO> ********** Densification in proj-0/systems/densification **********
INFO> Running Gromacs: minimization
INFO> Running Gromacs: nvt ensemble; 10.00 ps, 300.00 K
INFO> Running Gromacs: npt ensemble; 200.00 ps, 300.00 K, 10.00 bar
INFO> Current box side lengths: 3.127 nm x 3.127 nm x 3.127 nm
INFO> Density 829.27
INFO> Running-average-Density 726.35
INFO> Rolling-10-average-Density 828.52
INFO> Densified coordinates in proj-0/systems/densification/densified-npt.gro
Next comes a report of the precure:
INFO> **************** Precure in proj-0/systems/precure ****************
INFO> Running Gromacs: npt ensemble; 200.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 3.133 nm x 3.133 nm x 3.133 nm
INFO> Density 824.17
INFO> Running-average-Density 828.76
INFO> Rolling-10-average-Density 823.60
INFO> Annealing: 5 points for 2 cycles over 160 ps
INFO> Annealed coordinates in annealed.gro
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 3.137 nm x 3.137 nm x 3.137 nm
INFO> Density 820.71
INFO> Running-average-Density 830.07
INFO> Rolling-10-average-Density 829.90
Next we begin the CURE iterations:
INFO> ********* Connect-Update-Relax-Equilibrate (CURE) begins **********
INFO> Attempting to form 120 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 1 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 1 will generate 39 new bonds
INFO> Step "cure_relax" initiated on 39 distances (max 0.463 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.405 0.619
INFO> 2 0.338 0.555
INFO> 3 0.277 0.504
INFO> 4 0.230 0.454
INFO> 5 0.165 0.424
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 3.043 nm x 3.043 nm x 3.043 nm
INFO> Density 894.71
INFO> Running-average-Density 843.62
INFO> Rolling-10-average-Density 892.91
INFO> Iteration 1 current conversion 0.307 or 39 bonds
This first iteration shows that, with a search radius of 0.5 nm, HTPolyNet identified 39 allowable bonds. It then forms them and progresses through the relaxation stages until they are at their correct lengths. Finally it runs the post-iteration NPT MD equilibration, reporting the resulting box dimensions and density.
Next we proceed through CURE iterations 2 through 13:
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 2 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 2 will generate 34 new bonds
INFO> Step "cure_relax" initiated on 34 distances (max 0.493 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.445 0.660
INFO> 2 0.386 0.640
INFO> 3 0.326 0.551
INFO> 4 0.275 0.509
INFO> 5 0.214 0.458
INFO> 6 0.165 0.415
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 3.002 nm x 3.002 nm x 3.002 nm
INFO> Density 927.98
INFO> Running-average-Density 889.38
INFO> Rolling-10-average-Density 927.66
INFO> Iteration 2 current conversion 0.575 or 73 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 3 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 3 will generate 17 new bonds
INFO> Step "cure_relax" initiated on 17 distances (max 0.499 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.460 0.652
INFO> 2 0.407 0.631
INFO> 3 0.329 0.551
INFO> 4 0.269 0.511
INFO> 5 0.213 0.463
INFO> 6 0.167 0.411
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.975 nm x 2.975 nm x 2.975 nm
INFO> Density 950.93
INFO> Running-average-Density 917.79
INFO> Rolling-10-average-Density 944.20
INFO> Iteration 3 current conversion 0.709 or 90 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 4 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 4 will generate 11 new bonds
INFO> Step "cure_relax" initiated on 11 distances (max 0.465 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.413 0.622
INFO> 2 0.349 0.574
INFO> 3 0.277 0.512
INFO> 4 0.219 0.461
INFO> 5 0.156 0.408
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.995 nm x 2.995 nm x 2.995 nm
INFO> Density 930.82
INFO> Running-average-Density 922.40
INFO> Rolling-10-average-Density 939.98
INFO> Iteration 4 current conversion 0.795 or 101 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 5 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 5 will generate 6 new bonds
INFO> Step "cure_relax" initiated on 6 distances (max 0.495 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.440 0.660
INFO> 2 0.379 0.594
INFO> 3 0.333 0.573
INFO> 4 0.272 0.503
INFO> 5 0.212 0.454
INFO> 6 0.166 0.412
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.973 nm x 2.973 nm x 2.973 nm
INFO> Density 950.45
INFO> Running-average-Density 918.67
INFO> Rolling-10-average-Density 953.45
INFO> Iteration 5 current conversion 0.843 or 107 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 6 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 6 will generate 3 new bonds
INFO> Step "cure_relax" initiated on 3 distances (max 0.481 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.417 0.633
INFO> 2 0.369 0.590
INFO> 3 0.318 0.540
INFO> 4 0.271 0.498
INFO> 5 0.206 0.440
INFO> 6 0.162 0.399
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.969 nm x 2.969 nm x 2.969 nm
INFO> Density 954.32
INFO> Running-average-Density 930.29
INFO> Rolling-10-average-Density 951.54
INFO> Iteration 6 current conversion 0.866 or 110 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 7 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 7 will generate 3 new bonds
INFO> Step "cure_relax" initiated on 3 distances (max 0.472 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.411 0.639
INFO> 2 0.342 0.588
INFO> 3 0.282 0.522
INFO> 4 0.224 0.464
INFO> 5 0.161 0.414
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.988 nm x 2.988 nm x 2.988 nm
INFO> Density 935.63
INFO> Running-average-Density 892.22
INFO> Rolling-10-average-Density 941.46
INFO> Iteration 7 current conversion 0.890 or 113 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 8 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 8 will generate 1 new bond
INFO> Step "cure_relax" initiated on 1 distance (max 0.349 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.294 0.526
INFO> 2 0.251 0.470
INFO> 3 0.203 0.420
INFO> 4 0.157 0.400
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.977 nm x 2.977 nm x 2.977 nm
INFO> Density 945.79
INFO> Running-average-Density 922.02
INFO> Rolling-10-average-Density 942.38
INFO> Iteration 8 current conversion 0.898 or 114 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 9 begins ~~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 9 will generate 1 new bond
INFO> Step "cure_relax" initiated on 1 distance (max 0.459 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.393 0.627
INFO> 2 0.342 0.566
INFO> 3 0.281 0.522
INFO> 4 0.215 0.443
INFO> 5 0.160 0.402
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.990 nm x 2.990 nm x 2.990 nm
INFO> Density 933.71
INFO> Running-average-Density 902.46
INFO> Rolling-10-average-Density 924.82
INFO> Iteration 9 current conversion 0.906 or 115 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 10 begins ~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Radius increased to 0.75 nm
INFO> Iteration 10 will generate 2 new bonds
INFO> Step "cure_relax" initiated on 2 distances (max 0.578 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.506 0.734
INFO> 2 0.458 0.690
INFO> 3 0.401 0.634
INFO> 4 0.336 0.569
INFO> 5 0.273 0.503
INFO> 6 0.216 0.445
INFO> 7 0.163 0.414
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 3.008 nm x 3.008 nm x 3.008 nm
INFO> Density 916.97
INFO> Running-average-Density 897.35
INFO> Rolling-10-average-Density 928.33
INFO> Iteration 10 current conversion 0.921 or 117 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 11 begins ~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 11 will generate 1 new bond
INFO> Step "cure_relax" initiated on 1 distance (max 0.423 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.379 0.616
INFO> 2 0.314 0.531
INFO> 3 0.260 0.496
INFO> 4 0.199 0.437
INFO> 5 0.153 0.404
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.990 nm x 2.990 nm x 2.990 nm
INFO> Density 933.55
INFO> Running-average-Density 886.99
INFO> Rolling-10-average-Density 922.68
INFO> Iteration 11 current conversion 0.929 or 118 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 12 begins ~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Iteration 12 will generate 1 new bond
INFO> Step "cure_relax" initiated on 1 distance (max 0.469 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.408 0.646
INFO> 2 0.345 0.581
INFO> 3 0.275 0.503
INFO> 4 0.218 0.448
INFO> 5 0.157 0.409
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.979 nm x 2.979 nm x 2.979 nm
INFO> Density 943.83
INFO> Running-average-Density 923.67
INFO> Rolling-10-average-Density 939.60
INFO> Iteration 12 current conversion 0.937 or 119 bonds
INFO> ~~~~~~~~~~~~~~~~~~~~~~~ Iteration 13 begins ~~~~~~~~~~~~~~~~~~~~~~~
INFO> Bond search using radius 0.5 nm initiated
INFO> Radius increased to 0.75 nm
INFO> Iteration 13 will generate 1 new bond
INFO> Step "cure_relax" initiated on 1 distance (max 0.571 nm)
INFO> Stage Max-distance (nm) Max-1-4-distance (nm)
INFO> 1 0.505 0.711
INFO> 2 0.449 0.672
INFO> 3 0.394 0.622
INFO> 4 0.328 0.538
INFO> 5 0.278 0.526
INFO> 6 0.208 0.451
INFO> 7 0.155 0.393
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.976 nm x 2.976 nm x 2.976 nm
INFO> Density 946.62
INFO> Running-average-Density 920.90
INFO> Rolling-10-average-Density 937.37
INFO> Iteration 13 current conversion 0.945 or 120 bonds
This meets our desired cure of 95%, so now HTPolyNet proceeds to capping,
and not finding any cappable bonds, proceeds to the postcure:
INFO> ************************* Capping begins **************************
INFO> Capping will generate 0 new bonds
INFO> ********** Connect-Update-Relax-Equilibrate (CURE) ends ***********
INFO> *************** Postcure in proj-0/systems/postcure ***************
INFO> Annealing: 5 points for 2 cycles over 160 ps
INFO> Annealed coordinates in annealed.gro
INFO> Running Gromacs: npt ensemble; 100.00 ps, 300.00 K, 1.00 bar
INFO> Current box side lengths: 2.967 nm x 2.967 nm x 2.967 nm
INFO> Density 955.28
INFO> Running-average-Density 948.42
INFO> Rolling-10-average-Density 953.16
INFO> *********** Final data to proj-0/systems/final-results ************
INFO> ********************* HTPolyNet runtime ends **********************
This just tells us the final density and where the final results are found. Let’s now take a look at the final results.
A Note About FutureWarnings¶
Depending on your versions of numpy and pandas, your console may report messages like this:
/home/cfa/anaconda3/envs/mol-env/lib/python3.10/site-packages/numpy/core/fromnumeric.py:59: FutureWarning: 'DataFrame.swapaxes' is deprecated and will be removed in a future version. Please use 'DataFrame.transpose' instead.
return bound(*args, **kwds)
This happens because I use numpy.split() to partition the set of potential bonds into chunks that are processed in parallel through the bond-check algorithm, and numpy.split() calls the .swapaxes() method on its argument. Since the argument is a pandas DataFrame, and the pandas team is apparently unhappy about their swapaxes method being named swapaxes and prefers the name transpose (which I agree is a better name), they have the DataFrame.swapaxes method throw this little warning. It is nothing to worry about, as far as I can see. It will likely one day be resolved by the numpy team.