---
# System prepended metadata

title: 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
<!--- >  see local folder: /Users/lixinmeng/Desktop/working/md-scf/HyMD-systems-test/AzoTMA/AAMD --->

## step 2, Martini CG
>  follows the protocol from Martini website,[link](http://cgmartini.nl/index.php/tutorials-general-introduction-gmx5/parametrzining-new-molecule-gmx5)
> <!--- - own reference: $ cd /Volumes/xinmeng-lacie/phd_bar/main/cg-martini --->

### 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 



```























































