Difference between revisions of "Setting up (CHARMM)"

From Docswiki
Jump to navigation Jump to search
import>Csw34
import>Csw34
Line 122: Line 122:
 
'''AS OF 07/02/2012 THIS STEP DOES NOT WORK AS THE SCRIPT ONLY SUPPORTS AMBER - THIS WILL BE CORRECTED SOON! FOR NOW, JUST TAKE A COPY OF THE BELOW perm.allow FILE :)'''
 
'''AS OF 07/02/2012 THIS STEP DOES NOT WORK AS THE SCRIPT ONLY SUPPORTS AMBER - THIS WILL BE CORRECTED SOON! FOR NOW, JUST TAKE A COPY OF THE BELOW perm.allow FILE :)'''
   
You should get a perm.allow file like this:
+
You should get a ''perm.allow'' file like this:
 
<pre>
 
<pre>
 
4
 
4

Revision as of 17:28, 7 February 2012

Starting from a PDB

For this tutorial, we're going to assume you want to start from a PDB. This can raise some formatting issues, as the PDB format is rarely stuck to exactly, and even the slightest error can cause CHARMM to complain, and fail to read it in. If you have problems, you should compare your PDB file to the one we use here, and make the necessary changes yourself in vi. If you already have a .crd file to use, congratulations! You can safely skip the next section - although you will still need a PDB file to generate a perm.allow for your system automatically.

NOTE: if your system is NOT a standard peptide/protein, you will need to make the perm.allow file manually!

You can download the PDB we're using as the starting point for this tutorial here.

Converting a PDB into CHARMM .crd (CARD) format

In order to run GMIN, OPTIM or PATHSAMPLE using CHARMM, your input must be in CHARMM '.crd' or CARD format. You can do this using scripts form the MMTSB toolset, or just using CHARMM itself.

Here is some example CHARMM input to do this. You will need to change the set pardir line to point to where you have your copy of the svn repository. If you're working within CUC3, this will probably mean simply replacing my CRSID (csw34) with your own.

! CHARMM input to convert a PDB file into a CHARMM .crd (CARD) file
set pardir "/home/csw34/svn/CHARMM31/toppar"
set top "toph19_eef1_perm.inp"
set par "param19_eef1_perm.inp"

! BOMLev sets the level of warnings what do not cause the program to exit. -5 = very lax
BOMLev -5

! PRNLEV sets the ammount of output you get from CHARMM. 5 = high
PRNLEV +5

! Read standard topology and parameter files.
OPEN READ CARD UNIT 1 NAME @pardir/@top
READ RTF CARD UNIT 1
CLOSE UNIT 1

OPEN READ CARD UNIT 2 NAME @pardir/@par
READ PARAMETER CARD UNIT 2
CLOSE UNIT 2

! Generate the PSF for met-enk
READ SEQUence CARD
*
5
TYR GLY GLY PHE MET
GENErate FIRS NTER LAST CTER SETUp

! Read in the PDB file
OPEN UNIT 20 NAME metenk.pdb READ CARD
READ COOR UNIT 20 PDB
CLOSE UNIT 20

! Fill the internal coordinate table
IC FILL PRESERVE
IC PARAMETERS
IC PURGE
IC BUILD

! Write out the .crd (CARD) file
OPEN UNIT 20 NAME metenk.crd WRITE CARD
WRITE COOR UNIT 20 CARD
CLOSE UNIT 20

STOP

Assuming you have CHARMM in your path, and that you are in the directory containing metenk.pdb, you can run this simply using:

charmm < pdb2crd.inp

Hopefully this will run just fine, and you CHARMM output will end like this:

                   NORMAL TERMINATION BY NORMAL STOP
                    MAXIMUM STACK SPACE USED IS   72000
                    STACK CURRENTLY IN USE IS         0
                    MOST SEVERE WARNING WAS AT LEVEL  1
                    HEAP PRINTOUT-  HEAP SIZE  10240000
                    SPACE CURRENTLY IN USE IS         0
                    MAXIMUM SPACE USED IS          3866
                    FREE LIST
            PRINHP> ADDRESS:         1 LENGTH:  10240000 NEXT:         0

                    $$$$$ JOB ACCOUNTING INFORMATION $$$$$
                     ELAPSED TIME:     0.28  SECONDS 
                         CPU TIME:     0.00  SECONDS 

If you get a skull and crossbones, ask someone for help! If you wanted to do this for your own system, you would of course need to also modify the PSF section appropriately. You should now have a metenk.crd file, identical to this one.

That's it! We're pretty much good to go. As we're going to be using OPTIM later though, this is a good time to make a perm.allow file for your system...

Making a perm.allow file for later

The perm.allow file specifies which groups can be permuted when aligning endpoints of paths. For example, it ensures that all three H atoms of NH3 are treated as being the same, so NH3 rotation does not affect alignment, and two structures that differ only by a rotation of NH3 which permutes these atoms are correctly identified as being identical. To make a perm.allow file, we need a PDB that includes the termini being specified. This means we cannot use the original PDB without some modification. In this case, we need to change the N-terminal TYR residue to NTYR, and the C-terminal MET residue to CMET. We can easily do this in vi:

First, we make a copy of the original PDB to modify, and open it in vi:

cp metenk.pdb metenk_perm.pdb
vi metenk_perm.pdb

The fastest way to make the changes is using a 'find and replace' command. So, for the N-terminus:

:%s/TYR     1/NTYR    1/

and for the C-terminus:

:%s/MET     5/CMET    5/

You should see the changes have been made! We then save and exit:

:wq

ASIDE: we could also have done this using sed and regular expressions, like this:

sed 's/TYR     1/NTYR    1/' metenk.pdb | sed 's/MET     5/CMET    5/' > metenk_perm.pdb

Your metenk_perm.pdb file should be identical to this one. Now that we have the altered PDB, we can run a script to generate the perm.allow file for us:

~/svn/SCRIPTS/make_perm.allow/perm-pdb-ua.py metenk_perm.pdb CHARMM

AS OF 07/02/2012 THIS STEP DOES NOT WORK AS THE SCRIPT ONLY SUPPORTS AMBER - THIS WILL BE CORRECTED SOON! FOR NOW, JUST TAKE A COPY OF THE BELOW perm.allow FILE :)

You should get a perm.allow file like this:

4
3 0
 1 2 4
2 1
 8 10 9 11
2 1
 32 33 34 35
2 0
 47 48

The first line specifies that there are four permutable groups, in this case, NH3 of NTYR, atoms in the ring of NTYR, atoms in the ring of PHE and the COO- of CMET.

3 0
 1 2 4

Taking the first two as examples for the NH3 of NTYR, '3 0' says that there are three identical atoms, and that no other triples need move with them. In this case, these atoms are 1, 2 and 4.

2 1
 8 10 9 11

For the atoms in the ring of NTYR, two pairs of C atoms are identical (there is a mirror plane down the middle of the ring). If you were to flip one of these pairs however, you must also flip the other as the whole ring has been rotated. This is how we get to '2 1', with the 2 designating a pair, and the 1 that there is a second pair that must also be permuted at the same time. The four atoms below are the two pairs, 8 and 10, 9 and 11.