Assigning Point Charges to non-standard molecules for MD simulation.

For people doing classical MD often face the problem that the PDB file may contain a non-standard amino acid  / chromophore / ligand. The force field (CHARMM,OPLSAA,AMBER etc) provided by standard MD engines like Gromacs and NAMD can not tackle anything other than the handful of known amino-acids and solvents. Though to do a full blown computation and ascertain that the point charges are consistent might take a good while here I provide a fairly good way to parameterize ligands using the AMBER force field.

You will need to install (Antechamber)

A tutorial is provided in :

Example – Here I show the same steps for an organic molecule “triphenylamine”

1. Step 1 : Translate structure file to pdb/mol2 format.
It is key to understand that a pdb file only stores a set of xyz coordinates respect to a origin. A mol2 file also has the connections written in the file. Hence a mol2 file is more precise than a pdb or xyz file. In this case I have a ‘*.xyz’ file so I use “Openbabel” suite of softwares to convert it into mol2. It is important to note that one must manually check the output file to see if it has the connections written into it.

>babel [-i xyz] [-o mol2] babel.mol2
Files –, babel.mol2

2. Step 2 : Translate the coordinate to align the origin with the center of mass of the system.

>translate -i babel.mol2 -o babelcent.mol2 -f mol2 -c center -a1 0
Files – babelcent.mol2

3. Step 3 : Prepare Gaussian input file
>antechamber -fi mol2 -fo gcrt -i babelcent.mol2 -o babaelcent_gcrt.gau
Rename ‘babaelcent_gcrt.gau’ to ‘babaelcent_gcrt.gjf’

Make following modifications to header :
>#HF/6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) fopt iop(6/50=1)

File – babelcent_gcrt.gjf

4. Step 4 : Run Gaussian 03.
It is important to note that higher versions of Gaussian will not work with the present *.gjf file. PBS script will vary depending on the cluster. The one given below is a general example.

Estiamted time for 50 atoms is 9 hours with 12 procs.
Estiamted time for 40 atoms is 1 hour with 12 procs.

Submit with the following command on your local cluster
>qsub gaussian03.pbs
To check status of the run use ‘qstat’

Hint about the *.gjf file –

  • The “nproc=8” refers that the program me would like to use multiple of 8 processors.
  • The “mem=12GB” assigned during the run is unto 12GB.
  • It uses quantum computation at the level of Hatree Fock with 6-31G* basis sets and tight binding approximations. It also requires a full optimization. In case you are not sure about the structure is better to go with a full ‘fopt’ than a partial optimization ‘popt’.
  • To know more I found the file Gaussiantutorial.pdf pretty helpful.
  • Make sure that the file *.out ends with the line “Normal termination of Gaussian 03”
  • Since we will work with the “triphenylamine_g03.out” file hence forth it is better to make a copy of it “triphenylamine_g03_copy.out” and use that for further computations.

File – gaussian03.pbs, triphenylamine_g03.out

5. Step 5 : Takeout the ESP data from the “*.out” file.

>espgen -i triphenylamine_g03_copy.out -o tea_g03.esp
File – tea_g03.esp

6. Step 6 : Generate “*.ac” file.

In case you are using one conformation for your molecule only like in this example you can use the following method. The next step will generate a detail of bond, angle and dihedral in the AMBER format. The output file “” is what is most interesting for us.

>antechamber -fi gout -fo ac -i triphenylamine_g03_copy.out -o -c resp
>Check using : parmchk -i -f ac -o frcmod
File –

7. Step 7: Extract RESP Charges
Extract resp charge to a charge file tea_g03.crg from

>antechamber -fi ac -i -c wc -cf tea.crg
File – tea.crg

8.Step 8 : Combine *.crg with *.mol2 file

Read in tea.crg and babelcent.mol2 and generate an ac file

>antechamber -fi mol2 -i babelcent.mol2 -c rc -cf tea.crg -fo ac -o
File –

9. Step 9 : Reassign charges using a reference.

>prepgen -i -o tea_babel_resp_gaff.prepc -f car
File – tea_babel_resp_gaff.prepc, NEWPDB.PDB

For people running Gromacs

10) Install ‘‘. Details – README.txt.
11) Run ./ -i NEW.pdb
12) Compare file in NEWPDB.acpype to antechamber for validity.
13) Put charge from antechamber on acpype topology generated.

PDB Structure Factor Reader

It is a common practice to check if your MD computations have converged by matching the Experimental Structure factor with the structure factors provided in the PDB file. Each PDB file gives a value of dynamic structural factor in the 11th column of the PDB file, which can be matched to outputs of your

More on structure factor at – Structure Factor Wiki
Example – Fig 5 of the paper – Anisotropy of fluctuation dynamics of proteins with an elastic network model.

One can compute the same for there MD trajectories using g_rmsf in Gromacs suite of MD programmes.

Download File –
Step 1 : Rename pdb file of interest to ‘mypdb.pdb’
or, change the line in the perl code to your filename – $inputFile = ‘mypdb.pdb’;
Step 2 : Type – perl
Step 3 : xmgrace -nxy *.xvg (This will allow you to visulaise the *.xvg in xmgrace)