<p>
#!/usr/bin/python
# -*- coding: utf-8 -*-
#
import sys, os, subprocess, glob, linecache, string
 
# CHARMM executable
CHARMM = "/home/blub/.soft/usr/local/bin/charmm.xxlarge"
# Number of C Atom to start with
Cs = 2
# Number of last C Atom
Cl = 16
# Number of trajectory (DCD) files
N = 1
# Filename to store the final average S over all C-D
S = 'S.dat'
# trajectory files to read. Will be extended with "$N.dcd"
# Link all your dcd files into . and name the links membrane_npgt{1..N}.dcd
T = 'membrane_npgt'
# Directory to store results in
D = 'scd_proper'
# Edit the CHARMM config below to fit your needs
# Do not edit the "DIR" and "dcdname" var it is set via python var above
 
CHARMMcfg_form ='''
! the axial avg order parameter for an input C no.;
! N: Trajectory files -- C: chain carbon number @C
! http://www.charmm.org/ubbthreads/showflat.php?Cat=0&amp;Number=6581&amp;page=0&amp;vc=1
 
! dcdname:  membrane_npgt (will be extendet with "@K.dcd" --> 
 This e-mail address is being protected from spambots. You need JavaScript enabled to view it
 )
! segid:  SEGID of the Lipids
! nsteps:  number of Steps in the trajectory
! nlipids:  Number of Lipids
! DIR:    dir to store output in
! MinResid:  First resid of POPC (71)
! MaxResid:  Last resid of POPC (98)
! charmm N:40 C:2 < order_parameters.inp >&amp; scd/v2.out
 
set DIR $DT
set dcdname $TT
set nsteps 101000
set MinResidUP 15
set MaxResidUP 42
set MinResidLO 71
set MaxResidLO 98
 
calc nlipidsUP @MaxResidUP - @MinResidUP + 1
calc nlipidsLO @MaxResidLO - @MinResidLO + 1
calc nlipids @nlipidsUP + @nlipidsLO
 
set segid MEMB
 
bomblev -2
! Read in Topology and  Parameter files
 
open unit 1 card read name top_all27_prot_lipid.rtf
read RTF card unit 1
close unit 1
 
open unit 1 card read name par_all27_prot_lipid.prm
read PARA card unit 1
close unit 1
 
stream toppar_all27_lipid_cholesterol.str
 
!Read PSF &amp; COOR
open read card unit 10 name step5_assembly.psf
read psf  card unit 10
 
 
calc mxa @nlipids * 6 * 4
calc mxs 10 + @nlipids * 2 * 4
calc mxt @nsteps * @N
 
correl maxtim @MXT maxa @MXA maxs @MXS
enter sx zero
enter sy dupl sx
enter axy dupl sx
 
 
set resid @MinResidUP
set k 1
! TRAILING LETTERS IN THE 'H' ATOM NAMES REQUIRE @{C}, NOT @C
! ENTEr name VECT [X,Y,Z,R,XYZ] repeat(2x(atom-spec))
! atom-spec = segid  resid atom-name
label elpUP
! C3* == PALMITOYL
enter y@K vect z @segid @resid H@{C}X  @segid @resid C3@C
enter c@K vect r @segid @resid H@{C}X  @segid @resid C3@C
enter z@K vect z @segid @resid H@{C}Y  @segid @resid C3@C
enter d@K vect r @segid @resid H@{C}Y  @segid @resid C3@C
incr k by 1
incr resid by 1
if @k le @nlipidsUP goto elpUP
 
set resid @MinResidLO
! TRAILING LETTERS IN THE 'H' ATOM NAMES REQUIRE @{C}, NOT @C
! ENTEr name VECT [X,Y,Z,R,XYZ] repeat(2x(atom-spec))
! atom-spec = segid  resid atom-name
label elpLO
! C3* == PALMITOYL
enter y@K vect z @segid @resid H@{C}X  @segid @resid C3@C
enter c@K vect r @segid @resid H@{C}X  @segid @resid C3@C
enter z@K vect z @segid @resid H@{C}Y  @segid @resid C3@C
enter d@K vect r @segid @resid H@{C}Y  @segid @resid C3@C
incr k by 1
incr resid by 1
if @k le @nlipids goto elpLO
 
set k 1
label klp
calc u @K + 7
open unit @U file read name @{dcdname}@K.dcd
incr k by 1
if @k .le. @N goto klp
 
traj firstu 8 nunit @N
 
! LOOP OVER LIPIDS; NORM, CALC SCD, ACCUM
! S(CD) = (3*(z/r)**2 - 1)/2
! sr    = (1.5*(w@K / a@K)**2 - 0.5) * (1 / @nlipids)
! S=0,5(3cos**2TETA-1)
! cosTETA = Z/R = w@K/a@K
 
calc rLIPIDS 1. / @nlipids
 
set k 1
label clp
! C3* == PALMITOYL
mantim y@K ratio c@K
mantim y@K squa
mantim y@K mult 1.5
mantim y@K shift -0.5
mantim y@K mult @rLIPIDS
mantim sx add y@K
mantim z@K ratio d@K
mantim z@K squa
mantim z@K mult 1.5
mantim z@K shift -0.5
mantim z@K mult @rLIPIDS
mantim sy add z@K
incr k by 1
if @k le @nlipids goto clp
 
 
 
! Average both H
mantime axy add sx
mantime axy add sy
mantime axy mult 0.5
set averxy ?aver
 
edit sx veccod 3 delta 0.001 skip 1000 offset 1.
 
! MIXED CASE IS IGNORED; @C IS THE CHAIN CARBON NO.
open write unit 1 card name @{DIR}/
 This e-mail address is being protected from spambots. You need JavaScript enabled to view it
 
write sx unit 1 dumb time
* dumb
*
 
end
 
open write unit 1 card name @{DIR}/
 This e-mail address is being protected from spambots. You need JavaScript enabled to view it
 
write title unit 1
*PALMITOYL C:@C
*@averxy
*
stop
'''
CHARMMcfg_template=string.Template(CHARMMcfg_form)
# Define all substitutable strings in the template
CHARMMcfg = CHARMMcfg_template.substitute(DT=D, TT=T)
 
open("S.inp", "w").write(CHARMMcfg)
 
try:
 os.mkdir(D, 0755)
except OSError:
 print >>sys.stderr, "Dir exists, using existing dir."
 
runs = Cs
while runs <= Cl:
 print "Processing C%s" % (runs,)
 runCHARMM = subprocess.Popen([CHARMM, "N:%s" % (N,), "C:%s" % (runs,)], stdin=open("S.inp", "r"), stderr=open( "%s/C%s.err.out" % (D, runs,), "w"), stdout = open("%s/C%s.out" % (D, runs,), "w"))
 runCHARMM.communicate(CHARMMcfg)
 #remove empty err files
 err = open("%s/C%s.err.out" % (D, runs,)).read()
 if err == '':
 os.remove("%s/C%s.err.out" % (D, runs,))
 runs = runs + 1
 
os.remove("S.inp")
 
values = []
for file in glob.glob(D+'/a*'):
 values.append( float(linecache.getline(file, 2).strip()) )
 linecache.clearcache()
 
Sav = sum(values) / len(values)
print "Average Palmitoyl Order Parameter", Sav
open(D+"/"+S, 'w').write(str(Sav))
 
 

Attachments:

File	File size	Last Modified
order_parameters-popc.py	4 Kb	07/29/09 12:30