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) http://ambermd.org/antechamber/antechamber.html

A tutorial is provided in : http://ambermd.org/antechamber/example.html

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] triphenylamine.xyz [-o mol2] babel.mol2
Files – triphenylamine.xyz, 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 :
>%chk=molecule
>%nproc=8
>%mem=12GB
>#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 “tea_g03.ac” is what is most interesting for us.

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

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

>antechamber -fi ac -i tea_g03.ac -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 babel_resp.ac

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

9. Step 9 : Reassign charges using a reference.

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

**********************
For people running Gromacs

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