View source for If you need to change the number of atoms (e.g. making a united-atom charmm19 .crd file, or if atoms are missing)

Making a charmm19 .crd file from an experimental PDB structure is not so straightforward, I've found. The same can be true if the PDB file has missing atoms.  Apart from charmm bombing (!), symptoms of problems when you put a dodgy created .crd file into GMIN/OPTIM are messages about the internal coordinates being far from expected values in charmm sections, and RSQ being less than 10^-9 in EEF1 (if you're using that solvent model...)

Here are my notes on how to correct some common problems:

1. Make sure the atom and residue numbers start from 1 in the PDB file.  Some hacking of residue names may be required - see the charmm19 topology file to find the appropriate match if charmm can't recognise a name.

2. Sample section of the input for a charmm run:

<pre>
OPEN UNIT 20 FORM NAME protein_new.pdb READ
! read the sequence from the PDB file
READ SEQUence PDB UNIT 20
GENErate SETUP
IC PURGE

! read the coordinates from the PDB file
READ COOR PDB UNIT 20 SELE ALL END
CLOSE UNIT 20

! the command below puts in any missing hydrogens
HBUILD

! Fill the internal coordinate table from the PDB Cartesians
IC FILL

! For information, write out the current internals
OPEN UNIT 20 NAME afterfill.ic WRITE CARD
WRITE IC UNIT 20 CARD
CLOSE UNIT 20

! Fill in any missing internals with parameter values, then update the 
! Cartesians to reflect this.
IC PARA
COOR INIT
IC SEED 1 1 N 1 1 CA 1 1 C
IC BUILd

! write out both the current Cartesians and internal coordinates
OPEN UNIT 20 NAME coords.crd WRITE CARD
WRITE COOR UNIT 20 CARD
CLOSE UNIT 20
OPEN UNIT 20 NAME coords.ic WRITE CARD
WRITE IC UNIT 20 CARD
CLOSE UNIT 20

</pre>

3. Check that the capping O atoms at the C-terminus have not been placed on top of each other, when charmm uses standard internal coordinate values to place them.  If they have, then I have a small program, f90 code [[Media:euler.txt | here]], that rotates a specified point by a specified angle about an axis specified by two points, and outputs the resulting coordinates of the rotated point. Put just the x-y-z coordinates (no atom labels/numbers or anything) of the point to be rotated, then of any two points on the rotation axis (e.g C_alpha and C), and then the rotation angle in degrees, into a file named "points".  Run the program (./<executable name>) in the directory containing "points" and see what you get!  Excuse the lack of docu for this program ;-)  Then paste the new coordinates into your newly-made .crd file.  

After that, you should be about ready to go!

--[[User:jmc49|jmc49]] 17:23, 19 August 2009 (BST)