AzoTMA from AAMD to Martini & HyMD

# AzoTMA from AAMD to Martini & HyMD The key steps are: - [x] 1, prepare AAMD simulation results (Gromacs) - [x] 2, AAMD -> Martini CG MD - [ ] 3, Martini -> HyMD ## step 1, AAMD files - Thanks to Victoria, the AAMD files are already avaible  ## step 2, Martini CG > follows the protocol from Martini website,[link](http://cgmartini.nl/index.php/tutorials-general-introduction-gmx5/parametrzining-new-molecule-gmx5) >  ### 1. extract one molecule trajectory in AAMD ``` bash # inside the HyMD-systems-test $ cd AzoTMA # swith to working folder $ ls AAMD Martini $ cd AAMD $ gmx make_ndx -f AzoCIS_100ns_9nm.gro -o AzoMol1.ndx -> r 1 -> q $ gmx trjconv -f AzoCIS_100ns_9nm.xtc -s AzoCIS_100ns_9nm.tpr -o AzoCIS_Mol1.xtc -n AzoMol1.ndx -pbc mol -ur compact -center -b 0 -e 10000 -dt 10 -> 11 -> 11 $ gmx trjconv -f AzoCIS_Mol1.xtc -s AzoCIS_100ns_9nm.tpr -o AzoCIS_Mol1.gro -n AzoMol1.ndx -b 0 -e 0 -> 11 ###gmx editconf -f AzoCIS_Mol1.gro -o AzoCIS_Mol1.gro -box 5 5 5 -c # view single molecule trajectory $ vmd AzoCIS_Mol1.gro AzoCIS_Mol1.xtc ``` > repeat for the Trans ```bash $ echo 11 11 | gmx trjconv -f AzoTRANS_100ns_9nm.xtc -s AzoTRANS_100ns_9nm.tpr -o AzoTRANS_Mol1.xtc -n AzoMol1.ndx -pbc mol -ur compact -center -b 0 -e 10000 -dt 10 $ echo 11 | gmx trjconv -f AzoTRANS_Mol1.xtc -s AzoTRANS_100ns_9nm.tpr -o AzoTRANS_Mol1.gro -n AzoMol1.ndx -b 0 -e 0 $ vmd AzoTRANS_Mol1.gro AzoTRANS_Mol1.xtc ``` ### 2. implement A mapping via A mapping NDX file > Since we are dealing with with one small molecule, just construct the ndx file by hand: mapping-azo.ndx ```bash # continued from above $ vmd AzoTRANS_Mol1.gro > index 0 4 7 10 13 14 16 18 20 22 23 24 25 26 28 30 32 34 35 36 39 42 45 46 50 54 # group the atoms into different CG bead groups in mapping-azo-size2Ring-ScaleNgroupMass-size14.ndx $ cat mapping-azo-size2Ring-ScaleNgroupMass-size14.ndx [ 1 ] 1 5 [ 2 ] 8 11 [ 3 ] 14 17 [ 4 ] 15 19 [ 5 ] 21 23 [ 6 ] 24 25 [ 7 ] 26 29 [ 8 ] 27 31 [ 9] 33 35 [ 10 ] 36 37 [ 11 ] 40 43 [ 12 ] 46 47 [ 13 ] 46 51 [ 14 ] 46 55 ``` **Extra: sclae mass of N+ which is shared by three Groups** (may skip) ```bash # create demo itp, tpr $ cp diazo-cis.itp ORI_diazo-cis.itp $ cp diazo-trans.itp ORI_diazo-trans.itp # edit the diazo-cis.itp and diazo-trans.itp files # change the mass in the line: # 46 opls_288 0 AZT N 18 -0.080 14.00670 # from 14.00670 to 4.66890 e.g. 46 opls_288 0 AZT N 18 -0.080 4.66890 ; xinmeng scaled the mass by three # get the demo tpr with scaled mass; work with trans $ cp AzoTRANS_100ns_9nm.top AzoTRANS_Mol1.top $ vi AzoTRANS_Mol1.top ## ---> AZT 1 $ gmx grompp -v -f nvt.mdp -c AzoTRANS_Mol1.gro -p AzoTRANS_Mol1.top -o AzoTRANS_Mol1_SCALE.tpr -maxwarn 1000 ``` - the atom index in vmg starts from 0 and in gro/gmx starts from 1 - other mapping is possible, save for the future, if necessary ### 3. Generate mapped CG trajecoty from AA trajecotry ```bash $ seq 0 13 | gmx traj -f AzoTRANS_Mol1.gro -s AzoTRANS_Mol1_SCALE.tpr -oxt AzoTRANS_Mol1_CG.gro -n mapping-azo-size2Ring-ScaleNgroupMass-size14.ndx -com -ng 14 $ seq 0 13 | gmx traj -f AzoTRANS_Mol1.xtc -s AzoTRANS_Mol1_SCALE.tpr -oxt AzoTRANS_Mol1_CG.xtc -n mapping-azo-size2Ring-ScaleNgroupMass-size14.ndx -com -ng 14 # view output $ vmd AzoTRANS_Mol1_CG.gro AzoTRANS_Mol1_CG.xtc $ cp AzoTRANS_Mol1_CG.gro ../Martini/AzoTRANS.gro ``` > It is almost sure that we can feed the _gmx traj_ by the xtc and tpr of the whole assembly, instead of _Mol1_ system. The _ndx_ makes sure that only the atoms belong to the first MOL is read --- **Now trying to use an old-old protocol I used** ### 4. Prepare a basic topology of the AzoTrans molecule #### 4.1 use the brain, here we already have a draft at hand e.g. ![ ](https://i.imgur.com/CLRf0YJ.png) :::warning :warning: Later the SP0 is changed to SP1 ::: ``` bash # continued from above # switch to CG folder $ cd ../Martini # modify the atom name in the gro file $ vi AzoTRANS.gro # resname the atoname accordingly $ cat AzoTRANS.gro AzoTrans 14 1AZT CT1 1 4.816 4.442 5.420 1AZT CT2 2 4.803 4.484 5.162 1AZT CR1 3 4.747 4.441 4.937 1AZT CR2 4 4.684 4.614 4.839 1AZT CR3 5 4.660 4.431 4.742 1AZT N1 6 4.595 4.525 4.542 1AZT CR4 7 4.499 4.609 4.347 1AZT CR5 8 4.439 4.431 4.273 1AZT CR6 9 4.396 4.606 4.159 1AZT O1 10 4.271 4.504 3.996 1AZT CH 11 4.213 4.434 3.787 1AZT NH1 12 4.252 4.512 3.554 1AZT NH2 13 4.298 4.342 3.544 1AZT NH3 14 4.133 4.388 3.562 8.97629 8.97629 8.97629 ```` #### 4.2 Implement bond informationn using VMD and mol2 file > vmd AzoTRANS.gro ---> see the following youtube video & > get the file **AzoTRANS.mol2** {%youtube n8MEodmqces %} #### 4.3 Add charge information in mol2 file ```shell $ cp AzoTRANS.mol2 AzoTRANS_charge.mol2 ### edit accordingly $ cat AzoTRANS_charge.mol2 @<TRIPOS>MOLECULE AZT 14 15 1 0 0 SMALL 1 **** Energy = 0 @<TRIPOS>ATOM 1 CT1 48.1600 44.4200 54.2000 C 1 AZT 0.000000 2 CT2 48.0300 44.8400 51.6200 C 1 AZT 0.000000 3 CR1 47.4700 44.4100 49.3700 C 1 AZT 0.000000 4 CR2 46.8400 46.1400 48.3900 C 1 AZT 0.000000 5 CR3 46.6000 44.3100 47.4200 C 1 AZT 0.000000 6 N1 45.9500 45.2500 45.4200 N 1 AZT 0.000000 7 CR4 44.9900 46.0900 43.4700 C 1 AZT 0.000000 8 CR5 44.3900 44.3100 42.7300 C 1 AZT 0.000000 9 CR6 43.9600 46.0600 41.5900 C 1 AZT 0.000000 10 O1 42.7100 45.0400 39.9600 O 1 AZT 0.000000 11 CH 42.1300 44.3400 37.8700 C 1 AZT 0.100000 12 NH1 42.5200 45.1200 35.5400 N 1 AZT 0.300000 13 NH2 42.9800 43.4200 35.4400 N 1 AZT 0.300000 14 NH3 41.3300 43.8800 35.6200 N 1 AZT 0.300000 ``` ::: info One may try other topologies. What's done here is an initial try and mainly for illustration. :mega: ::: #### 4.3' Try to convert the mol2 to pdb file > My old protol use the Marterial Studio [commertial] to generate PDB with bonding info and mol2 file is only for charge; if MOL2 to PDB fails, then I will continue with MOL2 :smile: - open Babel should be a good option ```shell $ brew install open-babel ``` > There is some confliction with **hdf5-mpi** in my MAC. To avoid messing up my HyMD development environment, I managed to convert using the following website. http://chemdb.ics.uci.edu/cgibin/BabelWeb.py > The convertion is from **MOL2-Tripos Sybyl** to **PDB-Protein Databank PDB file** ```shell vmd AzoTRANS.pdb ``` - note in the coverted pdb file, still charge not included #### 4.4 obtain bonding ndx files :::info :fire: Starting now, it is a good practice to use our own parser to freely extract information. - following one old protool (I used in AAMD), I still use the MOL2 for the charge-info and PDB for bond-info. - Due to a 'trauma' I suffered about half year ago in an protocol of AAMD simulation, I define the atomnames in the MOL2/PDB with no repeat in the whole residue/molecule. Thus I also have to provide the atomtype information (see the last element in items from the handdrawing). This is down via replacing the *atomnames* by 'atomtypes' in the MOL2 file. ``` shell $ cat AzoTRANS_charge_atomtype.mol2 ``` - Here we try to get the ndx files; the ndx files will be used to get the reference bond/angle values/distribution. - Meanwhile, a 'dummy' itp is generated, where the bonding parametes are taken from the martini itp file. _Hopfully, later, the itp will be refined_ :monkey: ::: ```mermaid graph TD; A[AzoTRANS.pdb] --bond/atomname/atomid--> C{{Gen_ndx_itp.py}}; B[AzoTRANS_charge_atomtype.mol2] --atomtype--> C C --> D[AZT_bond.ndx]; C --> E[AZT_angle.ndx]; C --> H[AZT_dummy.itp] D --GMX--> F[(refence bond/angle )]; E --GMX--> F; ``` ```shell $ python Gen_ndx_itp.py ``` - The outcome flow chart inside the .py script (using our flow decorator) ![](https://i.imgur.com/Z5vZk9S.png) #### 4.5' Generate reference bond/angles ```shell # seq 0 14 | gmx distance -f ../AAMD/AzoTRANS_Mol1_CG.xtc -len 0.25 -tol 0.35 -n AZT_bond.ndx -oh AZT-bond-distribution.xvg -oall seq 0 14 | gmx distance -f ../AAMD/AzoTRANS_Mol1_CG.xtc -len 0.5 -tol 1 -n AZT_bond.ndx -oh AZT-bond-distribution.xvg -oall AZT-bond-data.xvg gmx angle -f ../AAMD/AzoTRANS_Mol1_CG.xtc -n AZT_angle.ndx -ov AZT-angle-data.xvg -type angle -all # ^-------- angle data file, first angle column is the average of all the angles ``` #### 4.6 Prepare for the dummy simulation > based on somefiles downloaded from http://cgmartini.nl/index.php/example-applications2/solvent-systems ```shell - prepare a 'big' solvent box gmx genconf -f nacl+wat.gro -o solv-box.gro -nbox 2 2 2 cat solv-box.gro --> box size 7.27924 7.27924 7.27924 - merge AZT inside gmx editconf -f AzoTRANS.gro -o AZT.gro -box 0 0 0 -c gmx insert-molecules -f solv-box.gro -ci AZT.gro -nmol 1 -o merge.gro -replace W -ip positions.dat -scale 1.5 - remove charge to make netral remove first two NA+ in the merge.gro file & change atomnum accordingly gmx editconf -f merge.gro -o merge.gro -resnr 1 @@ edit the box size line to be gapped with one space X X X - group residues together python simple-group-residues.py gmx editconf -f merge_ordered.gro -o merge_ordered.gro -resnr 1 > the sequence follows the 'tradition'; and easier for the ndx based analysis - edit topol.top cat topol.top ---------------- #include "martini_v2.2.itp" #include "martini_v2.0_ions.itp" #include "AZT_dummy.itp" [ system ] ; name AZT SOLV [ molecules ] ; name number AZT 1 W 3092 NA+ 62 CL- 64 ---------------- ``` #### 4.7 run > vi AZT_dummy.itp ===> SP0 --> SP1 gmx grompp -f minimization.mdp -c merge_ordered.gro -p topol.top -o em -maxwarn 100 gmx mdrun -v -deffnm em vmd em.gro gmx grompp -f relax.mdp -c em.gro -p topol.top -o relax -maxwarn 100 gmx mdrun -v -deffnm relax gmx grompp -f md.mdp -c relax.gro -p topol.top -o md -maxwarn 100 gmx mdrun -v -deffnm md :::warning In the martini cg, nrexcl = 1 in gromacs. May need to make it larger in dense/fine-grained systems.However rember in HyMD, there is no such exclusion operations. - maybe split the charges into a group of beads is too ***funny*** ... - but effectly, it is what we would like to achive, i.e. the bulky charge group!! ::: ---> #### 4.8 check bond (angle) distribution > in old simualtions, things are done in a py script; hard to pickup & restart now in Martini folder: seq 0 14 | gmx distance -f md.xtc -len 0.5 -tol 1 -n AZT_bond.ndx -oh AZT-bond-distribution-cg.xvg -oall AZT-bond-data-cg.xvg > in our script, source distance data is needed, thus AZT-bond-data-cg.xvg is gmx angle -f md.xtc -n AZT_angle.ndx -ov AZT-angle-data-cg.xvg -type angle -all python simple-plot-distribution.py - the dummy-simulation bond-distribution plots {%youtube fyhsLA0Ikvg %} - Modify by hand of the itp file cp AZT_dummy.itp AZT_dummy2.itp cp topol.top topol2.top ; change the new itp bash run.sh ; include the MD run and plot !!! Usually we could change the bond equilibrium distance or bond strength to have a better fitting. !!! Here, however, we are using a quite re-fined model. The shift is mainly due to the neibouring non-bonded exclusion interactions. As by default the nrexcl in the molecule itp is set to 1 !!! This neigbouring exludiosn could be a problem the non-bonded potential is hard. A solution is just change the nrexcl to 2. In hymd model, the field potential is soft, thus the bonded terms should be ok, even there is no nrexcl options. A. set nrexcl=2 B. further adjust the cg equlibrium according to bond distribution values ===> see the result: {%youtube 6iaFlSlu67M %} ====== angle distribution ====== edit according to the aa-distribution data see the result: {%youtube J8hUPxI3-c0 %} #### 4.9 prepare for a SA simulation > follows the setup in the AzoTRANS_100ns_9nm.top - the simulation consists AZT 90 SOL 22313 BR 90 - requred files: - AZT.gro ; AZT molecule - to be inserted in water box, SA target prepare 90: gmx insert-molecules -ci AZT.gro -o AZT90.gro -nmol 90 -box 9 9 9 -seed 1 - topolSA90.top ; top file, record the correct number of compoment molecules cp topolSA90.top topolSA.top - water.gro ; small equilibrated water box in martini; to fill in the big box solvate the AZT90.gro with water: gmx solvate -cp AZT90.gro -cs water.gro -o with_sol.gro -p topolSA.top -scale 1 - ION ;;;- cl.gro ; CL- ion; aa is Br ion; centered at 0 ; - to be inserted in water box to nutralize the system ;;; CAN NOT READ TOP so not using insert-moleucles now (gmx insert-molecules -f with_sol.gro -ci cl.gro -nmol 90 -o sat0.gro -replace W -scale 1 -p topolSA.top) cp minimization.mdp ions.mdp gmx grompp -f ions.mdp -c with_sol.gro -p topolSA.top -o ions.tpr -maxwarn 100 gmx genion -s ions.tpr -o allSA.gro -p topolSA.top -pname NA+ -nname CL- -neutral - local run - single case test gmx grompp -f minimization.mdp -c allSA.gro -p topolSA.top -o emSA -maxwarn 100 gmx mdrun -v -deffnm emSA vmd emSA.gro gmx grompp -f relax.mdp -c emSA.gro -p topolSA.top -o relaxSA -maxwarn 100 gmx mdrun -v -deffnm relaxSA gmx grompp -f md001dt.mdp -c relaxSA.gro -p topolSA.top -o mdSA -maxwarn 100 gmx mdrun -v -deffnm mdSA vmd relaxSA.gro mdSA.xtc ### 5. Run a set of SA simulations > on saga or betzy bash gulugulu-mix.sh ## step 3, HyMD > for the simplicity, generate the GMX thing by hand .... :::warning - The $\chi$ information has to be given e.g. in a forcefield.itp file - The bond/angle parameters, unfortunately, are not explicitely treated via pairs, but by their bead types ... This means that bond/angle parameters in the molecule itp can not be read/used, which is NOT NICE !!! Maybe possible to fixxxxxxxx this, e.g. in the force.py. For now maybeeeee, extract what is needed to able to have a gooooooooo ::: - Give the $\chi$ information > - not using the ML optimization, purily use the simple F-H model > - editing an martini_hymd_forcefield.itp (copied from Martini) ```bash cd hPF cp ../Martini/AZT_dummy2.itp . ; clean the itp file ```