CHARMM Development Project Forums CHARMM Community Script Archive lipid C-D order parameters using CORREL

The following 3 scripts compute the carbon-deuterium order parameter, S(C-D), for all of the aliphatic C-H bond vectors in the lipid DPPC. The main assumption is that the lipids are oriented, i.e. in a bilayer with the normal aligned with the Z axis (a pre-requisite for bilayer simulations). An assumed consequence is that lipids in a bilayer rotate around their respective long axis easily, but rarely invert, which would usually involve switching leaflets. The X and Y axis projections of bond vectors are expected to be averaged due to the rotation around the long axis, and the formula for the order parameter simplifies to

S(CD) = (3*z**2 - 1)/2

where z is the Z axis projection of the normalized bond vector (a unit vector). Much of the complexity of the script is bookkeeping related to the fact that DPPC has 80 C-H bond vectors to consider, and that one must consider all of the lipid molecules in the simulation system. I broke it down into 3 scripts to keep each simpler, and for efficiency reasons as well; to compute S(CD) for chain position 5 for an 80 lipid system requires applying the above equation to 320 vectors.

Nvect = (2 C-H bonds) * (2 chains) * (80 lipids) = 320

One script computes S(CD) for the CH2 groups of the alpha chain (the ethanolamine fragment), a second does the 5 vectors for the glycerol, and and the third computes it for the 4 C-H vectors at a specific chain position. The chain script requires an extra command line argument, which is the the chain position number (2-15 for DPPC); all the scripts expect an argument for the number of files. Assuming 40 .trj files (also called .dcd files) and input filenames of scd-eth.inp, scd-gly.inp, and scd-chn.inp

charmm N:40 < scd-eth.inp >& scd-eth.out
charmm N:40 C:5 < scd-chn.inp >& scd-chn5.out

scd-eth.inp; extra comments

* axial avg order parameter for ethanolamine C:H vectors; @N files
*

bomblev -1

! READ STD FILES FOR PROJECT; RTF, PARAM, PSF, COOR
stream rtfprm.str
stream psfcrd.str

! COMPUTE CORREL LIMITS
calc mxa 80 * 6 * 4
calc mxs 10 + 80*2 * 4
calc mxt 200 * @N

! INVOKE CORREL
correl maxtim @MXT maxa @MXA maxs @MXS
! ACCUMULATE AVERAGES IN THESE TIME SERIES; DECLARE IN THIS ORDER
enter sr zero
enter ss dupl sr
enter sx dupl sr
enter sy dupl sr
! LOOP OVER THE LIPIDS; USE ONLY Z AND R OF VECTOR
set k 1
label elp
enter w@K vect z L @K H11A  L @K C11
enter a@K vect r L @K H11A  L @K C11
enter x@K vect z L @K H11B  L @K C11
enter b@K vect r L @K H11B  L @K C11
enter y@K vect z L @K H12A  L @K C12
enter c@K vect r L @K H12A  L @K C12
enter z@K vect z L @K H12B  L @K C12
enter d@K vect r L @K H12B  L @K C12
incr k by 1
if k le 80 goto elp

! START FROM UNIT 8; ALLOWS UP TO CA. 90 FILES
set k 1
label klp
calc u @K + 7
open unit @U file read name dyn@K.trj
incr k by 1
if k .le. @N goto klp
traj firstu 8 nunit @N

! LOOP OVER 80 LIPIDS; NORM, CALC SCD, ACCUM
calc r80 1. / 80.
set k 1
label clp
mantim w@K ratio a@K   ! DIVIDE VECTOR Z BY R (z OF UNIT VECTOR)
mantim w@K squa        ! z**2
mantim w@K mult 1.5    ! * 3/2
mantim w@K shift -0.5  ! - 1/2
mantim w@K mult @R80   ! * 1/Nlipid
mantim sr add w@K      ! ACCUMLATE, H11A
mantim x@K ratio b@K
mantim x@K squa
mantim x@K mult 1.5
mantim x@K shift -0.5
mantim x@K mult @R80
mantim ss add x@K      ! ACCUMLATE, H11B
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 @R80
mantim sx add y@K      ! ACCUMLATE, H12A
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 @R80
mantim sy add z@K      ! ACCUMLATE, H12B
incr k by 1
if k le 80 goto clp
! LOOP OVER LIPIDS ENDS; COMBINE THE FOUR TIME SERIES
edit sr veccod 4 delta 0.001 skip 1000 offset 1.
! STASH DATA IN A SUBDIR NAMED scd; 5 COLS, TIME IN 1st COL
open write unit 1 card name scd/eth.dat
write sr unit 1 dumb time
* dumb
*
end
stop

scd-gly.inp; see comments in scd-eth.inp

* the axial avg order parameter for glycerol C:H vectors; @N files
*

bomblev -2
stream rtfprm.str
stream psfcrd.str

calc mxa 80 * 6 * 5
calc mxs 10 + 80*2 * 5
calc mxt 200 * @N

correl maxtim @MXT maxa @MXA maxs @MXS
enter sa zero
enter sb dupl sa
enter ss dupl sa
enter sx dupl sa
enter sy dupl sa
set k 1
label elp
enter v@K vect z L @K HS  L @K C2
enter r@K vect r L @K HS  L @K C2
enter w@K vect z L @K HA  L @K C1
enter a@K vect r L @K HA  L @K C1
enter x@K vect z L @K HB  L @K C1
enter b@K vect r L @K HB  L @K C1
enter y@K vect z L @K HX  L @K C3
enter c@K vect r L @K HX  L @K C3
enter z@K vect z L @K HY  L @K C3
enter d@K vect r L @K HY  L @K C3
incr k by 1
if k le 80 goto elp

set k 1
label klp
calc u @K + 7
open unit @U file read name dyn@K.trj
incr k by 1
if k .le. @N goto klp
traj firstu 8 nunit @N
! LOOP OVER 80 LIPIDS; NORM, CALC SCD, ACCUM
calc r80 1. / 80.
set k 1
label clp
mantim w@K ratio a@K
mantim w@K squa
mantim w@K mult 1.5
mantim w@K shift -0.5
mantim w@K mult @R80
mantim sa add w@K
mantim x@K ratio b@K
mantim x@K squa
mantim x@K mult 1.5
mantim x@K shift -0.5
mantim x@K mult @R80
mantim sb add x@K
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 @R80
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 @R80
mantim sy add z@K
mantim v@K ratio r@K
mantim v@K squa
mantim v@K mult 1.5
mantim v@K shift -0.5
mantim v@K mult @R80
mantim ss add v@K
incr k by 1
if k le 80 goto clp
edit sa veccod 5 delta 0.001 skip 1000 offset 1.
open write unit 1 card name scd/gly.dat
write sa unit 1 dumb time
* dumb
*
end

scd-chn.inp; bonus comments included; see also scd-eth.inp

* the axial avg order parameter for an input C no.; @N files
* additional required arg C: chain carbon number @C
*

bomblev -2
stream rtfprm.str
stream psfcrd.str

calc mxa 80 * 6 * 4
calc mxs 10 + 80*2 * 4
calc mxt 200 * @N

correl maxtim @MXT maxa @MXA maxs @MXS
enter sr zero 
enter ss dupl sr
enter sx dupl sr
enter sy dupl sr
set k 1
! TRAILING LETTERS IN THE 'H' ATOM NAMES REQUIRE @{C}, NOT @C
label elp
enter w@K vect z L @K H@{C}R  L @K C2@C
enter a@K vect r L @K H@{C}R  L @K C2@C
enter x@K vect z L @K H@{C}S  L @K C2@C
enter b@K vect r L @K H@{C}S  L @K C2@C
enter y@K vect z L @K H@{C}X  L @K C3@C
enter c@K vect r L @K H@{C}X  L @K C3@C
enter z@K vect z L @K H@{C}Y  L @K C3@C
enter d@K vect r L @K H@{C}Y  L @K C3@C
incr k by 1
if k le 80 goto elp

set k 1
label klp
calc u @K + 7
open unit @U file read name dyn@K.trj
incr k by 1
if k .le. @N goto klp
traj firstu 8 nunit @N
! LOOP OVER LIPIDS; NORM, CALC SCD, ACCUM
calc r80 1. / 80.
set k 1
label clp
mantim w@K ratio a@K
mantim w@K squa
mantim w@K mult 1.5
mantim w@K shift -0.5
mantim w@K mult @R80
mantim sr add w@K
mantim x@K ratio b@K
mantim x@K squa
mantim x@K mult 1.5
mantim x@K shift -0.5
mantim x@K mult @R80
mantim ss add x@K
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 @R80
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 @R80
mantim sy add z@K
incr k by 1
if k le 80 goto clp
edit sr veccod 4 delta 0.001 skip 1000 offset 1.
! MIXED CASE IS IGNORED; @C IS THE CHAIN CARBON NO.
open write unit 1 card name scd/v@C.dat
write sr unit 1 dumb time
* dumb
*
end
stop

Dear Dr Venable,

Please correct my understanding.

In case of plotting between scd and position of the carbon atoms, I need to make average; say sr and ss for C11 and so on.

Is it correct?

Thank you.

In general, yes, the Scd for each of the 2 C-H vectors on a given C atom are averaged, as they cannot be distinguished experimentally in most cases. The exception is C2 of each chain (C22 and C32), especially the beta chain (C22); the C-H vectors are inequivalent and can be distinguished by the NMR quadrapole splitting experiment. Computing the Scd for the 2 vectors separately is a test of convergence; if <Scd> (average over time) for 2 C-H vectors at a given chain position don't agree that well, the simulation may need to be run longer.

Let me ask you a question on S(CD) parameter.
In order to compare to experimental results, the simulation setup should be exactly the same as in experiments, e.g, the same temperature, pressure, surface tension ?

Thanks.

Certainly, the temperature, hydration level, and the surface area per lipid need to be fairly close, as the order parameters are known to depend on these properties.

Which one is the most important factor among the above three ?
Is it necessary to check the convergence of the order parameter ?

Thanks.

Dear Prof. Rick,

I have to analyse the order parameter for more than 100 files, the script mentioned here can take maximum of 90 files. Is there any way of reading more than 100 trj files. I have tried the close unit @u also but that too doesnt seem to be working.

waiting for your reply,
regards,
MzZt.

Update: CHARMM now supports many open files, whatever the OS will allow; start from UNIT 101 and increment upwards.

You have to use the MERGE command (dynamc.doc) to combine files; I often do a temporary merge from the storage location to the /scratch volume of a machine, then run a series of analysis on files that have been merged to 1 ns chunks. Since this is the Script Archive, here's the script; the @J variable must be passed from the command line--

* merge 10:1
*

bomblev -1

! READ RTF, PARAM FILES
stream rtfprm.str
! READ PSF, COOR
stream psfcrd.str
! AND SETUP CRYSTAL (IMPORTANT!!!)
stream cryst.str

! COMPUTE LAST OF SET
calc ilst @J + 9
! START FROM UNIT 21
set iu 21
! INPUT FILE OPEN LOOP
label flp
open unit @IU read file name dyn@J.trj
incr iu by 1
incr j by 1
if j .le. @ILST goto flp

! OPEN THE OUTPUT FILE AS UNIT 19
open unit 19 write file name /scratch/dyn@K.trj
! MERGE THE FILES TO CREATE A NEW ONE
merge coor firstu 21 nunit 10 output 19 nfile 1000

stop

An added script that illustrates detection of C=C double bonds via the CEL1 chemical atom type; 72 DOPC lipids, segid D.

* the axial avg order parameter for an input C no. @C
* modified for C=C case; first file @FF last file @LF
* subdir @S
*

stream rtfprm.str
stream psfcrd.str

! number of lipids
set nl = 72
! number of files
calc n = 1 + ( @LF - @FF )

calc mxa @NL * 6 * 4
calc mxs 10 + @NL*2 * 4
calc mxt 500 * @N

! unsaturation check; beta chain
set ub = 0
define unsat sele atom D 1 C2@C .and. chem CEL1 end
if ub lt ?NSEL set ub ?NSEL
! unsaturation check; gamma chain
set ug = 0
define unsat sele atom D 1 C3@C .and. chem CEL1 end
if ug lt ?NSEL set ug ?NSEL

correl maxtim @MXT maxa @MXA maxs @MXS
enter sr zero
if ub eq 0 enter ss dupl sr
enter sx dupl sr
if ug eq 0 enter sy dupl sr
set k 1
label elp
enter w@K vect z D @K H@{C}R D @K C2@C
enter a@K vect r D @K H@{C}R D @K C2@C
if ub eq 0 enter x@K vect z D @K H@{C}S D @K C2@C
if ub eq 0 enter b@K vect r D @K H@{C}S D @K C2@C
enter y@K vect z D @K H@{C}X D @K C3@C
enter c@K vect r D @K H@{C}X D @K C3@C
if ug eq 0 enter z@K vect z D @K H@{C}Y D @K C3@C
if ug eq 0 enter d@K vect r D @K H@{C}Y D @K C3@C
incr k by 1
if k le @NL goto elp

! OPEN FILES, FILL SERIES VIA TRAJ
set k = @FF
set u = 8
label opnlp
open read unit @U file name @S/dyn@K.trj
incr u by 1
incr k by 1
if k le @LF goto opnlp
traj firstu 8 nunit @N

! LOOP OVER LIPIDS; NORM, CALC SCD, ACCUM
set k 1
label clp
mantim w@K ratio a@K
mantim w@K squa
mantim w@K mult 1.5
mantim w@K shift -0.5
mantim sr add w@K
if ub eq 0 then
mantim x@K ratio b@K
mantim x@K squa
mantim x@K mult 1.5
mantim x@K shift -0.5
mantim ss add x@K
endif
mantim y@K ratio c@K
mantim y@K squa
mantim y@K mult 1.5
mantim y@K shift -0.5
mantim sx add y@K
if ug eq 0 then
mantim z@K ratio d@K
mantim z@K squa
mantim z@K mult 1.5
mantim z@K shift -0.5
mantim sy add z@K
endif
incr k by 1
if k le @NL goto clp

! use final mantime to set AVER, FLUC for output
open write unit 1 card name scd/chn@{C}scd.txt
echu 1
mantim sr divi @NL
echo sr@C ?AVER ?FLUC
if ub eq 0 then
mantim ss divi @NL
echo ss@C ?AVER ?FLUC
endif
mantim sx divi @NL
echo sx@C ?AVER ?FLUC
if ug eq 0 then
mantim sy divi @NL
echo sy@C ?AVER ?FLUC
endif
close unit 1

! *** uncomment the next lines to write out the time series
! edit sr veccod 4 delta 0.001 skip 1000 offset 1.
! open write unit 1 card name scd/vchn@C.dat
! write sa unit 1 dumb time
! * dumb
! *

end
stop

A similar script to the one above for C=C in DOPC, but modified for POPC. For some reason, the POPC vinyl H atoms break the standard naming conventions for the oleate chain. Also, an exception was added to allow for the longer beta chain.

The first and last parts of script have been omitted; see the DOPC example above for the missing bits.

correl maxtim @MXT maxa @MXA maxs @MXS
enter sr zero
if ub eq 0 enter ss dupl sr
enter sx dupl sr
if ug eq 0 enter sy dupl sr
set k 1
label elp
if ub eq 0 then
enter w@K vect z D @K H@{C}R D @K C2@C
enter a@K vect r D @K H@{C}R D @K C2@C
enter x@K vect z D @K H@{C}S D @K C2@C
enter b@K vect r D @K H@{C}S D @K C2@C
else
enter w@K vect z D @K H@{C}1 D @K C2@C
enter a@K vect r D @K H@{C}1 D @K C2@C
endif
if c gt 15 goto eskip
if ug eq 0 then
enter y@K vect z D @K H@{C}X D @K C3@C
enter c@K vect r D @K H@{C}X D @K C3@C
enter z@K vect z D @K H@{C}Y D @K C3@C
enter d@K vect r D @K H@{C}Y D @K C3@C
else
enter y@K vect z D @K H@{C}1 D @K C3@C
enter c@K vect r D @K H@{C}1 D @K C3@C
endif
label eskip
incr k by 1
if k le @NL goto elp

! OPEN FILES, FILL SERIES VIA TRAJ
set k = @FF
set u = 101
label opnlp
open read unit @U file name @S/dyn@K.trj
incr u by 1
incr k by 1
if k le @LF goto opnlp
traj firstu 101 nunit @N

! LOOP OVER LIPIDS; NORM, CALC SCD, ACCUM
set k 1
label clp
mantim w@K ratio a@K
mantim w@K squa
mantim w@K mult 1.5
mantim w@K shift -0.5
mantim sr add w@K
if ub eq 0 then
mantim x@K ratio b@K
mantim x@K squa
mantim x@K mult 1.5
mantim x@K shift -0.5
mantim ss add x@K
endif
if c gt 15 goto mskip
mantim y@K ratio c@K
mantim y@K squa
mantim y@K mult 1.5
mantim y@K shift -0.5
mantim sx add y@K
if ug eq 0 then
mantim z@K ratio d@K
mantim z@K squa
mantim z@K mult 1.5
mantim z@K shift -0.5
mantim sy add z@K
endif
label mskip
incr k by 1
if k le @NL goto clp

! use final mantime to set AVER, FLUC for output
open write unit 1 card name chn@{C}scd@FFI-@LFI.txt
echu 1
mantim sr divi @NL
echo sr@C ?AVER ?FLUC
if ub eq 0 then
mantim ss divi @NL
echo ss@C ?AVER ?FLUC
endif
if c gt 15 goto oskip
mantim sx divi @NL
echo sx@C ?AVER ?FLUC
if ug eq 0 then
mantim sy divi @NL
echo sy@C ?AVER ?FLUC
endif
label oskip
close unit 1

Hello,

What does "Error in parameter substitution" mean and how can it be fixed?

Thank you,
Mihaela

It probably means you've made a mistake.

Unless your question is explicitly about the scripts posted above, this may also be the wrong forum for this post. I suggest you review the guidelines in the READ BEFORE POSTING topic more carefully.

The question was about POPC script above (page 2). Could you please tell me how to define chain carbon number @C?

CHARMM> define unsat sele atom D 1 C2@C .and. chem CEL1 end
* WARNING * Command ignored. Token not found: >C<

***** LEVEL 0 WARNING FROM *****
***** Error in parameter substitution.

Thank you,
Mihaela

The value for @C is passed as a command line argument; see the very first post in this thread for an example.

Yes, I do that, and the next error I get is:

ERROR IN NXTATM: Unrecognizable SEGID or residue number "D ".

This error comes after reading rtf, prm, psf and crd files. How do I identify what D is?
Thank you,
Mihaela

As indicated in the error message, D is expected to be a SEGID name; that was the name I used for the lipid bilayer segment in my simulation. You must change that to the name used in your bilayer model. Most of the scripts here in the Archive will require some editing to work properly.

See also select.doc, which describe CHARMM's powerful atom selection feature.

I found it, thank you. Now i don't get errors in output file. Only after I run the program, I get:

At line 145 of file /home/themis/c39a2/source/util/parse.src (unit = 5, file = 'stdin')
Fortran runtime error: Sequential READ or WRITE not allowed after EOF marker, possibly use REWIND or BACKSPACE

What does it mean?

I'm guessing that there's a problem with the input and some clues about the problem in the output; without more information it's difficult to make any suggestions.

Also, could you please tell me if these numbers:

calc mxa @NL * 6 * 4
calc mxs 10 + @NL*2 * 4
calc mxt 500 * @N

will have to change in my script?

What is the significance if 15 here:

if c gt 15 goto eskip
if ug eq 0 then

Thank you,
Mihaela

Also, in the line:

define unsat sele atom D 1 C2@C .and. chem CEL1 end

what is the meaning of 1, chem and CEL1?
Thank you,
Mihaela

calc mxa @NL * 6 * 4
calc mxs 10 + @NL*2 * 4
calc mxt 500 * @N

These lines compute the values needed for CORREL startup, i.e. MAXAtoms, MAXSeries, and MAXTimesteps; their meaning is described in correl.doc

The value '500' may need to be changed, as it is based on the number of stored coordinate sets in each trajectory file. The variable @NL is the number of lipid molecules, and is set near the top of the main script.


if c gt 15 goto eskip
if ug eq 0 then

The above two lines are exception code added specifically for POPC, where the gamma chain (a.k.a. chain 1) has only 16 C atoms.


define unsat sele atom D 1 C2@C .and. chem CEL1 end
what is the meaning of 1, chem and CEL1?

Another atom selection issue (select.doc); the CHEM token allows one to select atoms based on their parameter type, rather than a IUPAC style label. The CEL1 parameter type is specific for a C atom involved in a C=C double bond.

When I run the program, I get no errors, but no atoms are selected:

CHARMM> ! unsaturation check; beta chain
CHARMM> set ub = 0
Parameter: UB <- "0"

CHARMM> define unsat sele atom A 1 C2@C .and. chem CEL1 end
Parameter: C -> "2"
SELRPN> 0 atoms have been selected out of 36716

What am I doing wrong?
Thank you,
Mihaela

Nothing is wrong; that test only succeeds if chain position 2 has a double bond (chem CEL1). For POPC, only positions 9 and 10 of the beta chain (chain 2) have a double bond.

Right, I get no errors in the output file, but the order parameter output is invariably 0, regardless what value I assign to C:

For example, for C=2, I get:

SR2 0 0
SS2 0 0
SX2 0 0
SY2 0 0

Where could the error be?

There are probably messages in the output file which give some indication of the problem.

I think the problem is in the trajectory (dcd file)for POPC. I had POPC in water, so I needed only the trajectory for POPC. Maybe my method of selecting only the lipids was wrong, because I got a mismatch in the number of atoms from psf file for lipids and the trajectory file for lipids.
Do you have a recommendation of how to obtain the trajectory file for POPC only from the trajectory file of POPC in water?

I usually keep the water.

If I keep the water, I get the same result:
SR2 0 0
SS2 0 0
SX2 0 0
SY2 0 0

What is to be done?

Find the clues in your output file.

OK, now for the final part, in the output file I get:

ORREL> ! use final mantime to set AVER, FLUC for output
CORREL> open write unit 1 card name cd@{C}.txt
Parameter: C -> "2"
At line 145 of file /home/themis/c39a2/source/util/parse.src (unit = 5, file = 'test.inp')
Fortran runtime error: Sequential READ or WRITE not allowed after EOF marker, possibly use REWIND or BACKSPACE
OPNLGU> Unit already open. The old file will be closed first.
VCLOSE: Closing unit 1 with status "KEEP"
VOPEN> Attempting to open::cd2.txt::
OPNLGU> Unit 1 opened for WRITE access to cd2.txt

I think this is the clue, but how do I interpret it?

The Fortran runtime error message suggests some kind of problem with the file named test.inp; perhaps that file has become corrupted in some way. You have not provided enough information to say any more.

Here is the file test.inp, for 274 POPC lipids.(I made it following the example above, with a saved trajectory file for lipids without water and a psf file written in charm format.The number of atoms in trajectory file match those in psf ):

! the axial avg order parameter for an input C no. @C
! modified for C=C case; first file @FF last file @LF
! subdir @S


! Read initial data files
open unit 1 read form name top_all36_lipid.rtf
read rtf card unit 1
close unit 1
open unit 1 read form name par_all36_lipid.prm
read param card unit 1
close unit 1

! Read parameter file in command line

open read card unit 1 name popc.psf
read psf card unit 1
open read card unit 1 name popc.crd
read coor card unit 1

! number of lipids
set nl = 274
! number of files
calc n = 1

calc mxa @NL * 6 * 4
calc mxs 10 + @NL*2 * 4
calc mxt 2000 * @N

! unsaturation check; beta chain
set ub = 0
define unsat sele atom A 1 C2@C .and. chem CEL1 end
if ub lt ?NSEL set ub ?NSEL
! unsaturation check; gamma chain
set ug = 0
define unsat sele atom A 1 C3@C .and. chem CEL1 end
if ug lt ?NSEL set ug ?NSEL

correl maxtim @MXT maxa @MXA maxs @MXS
enter sr zero
if ub eq 0 enter ss dupl sr
enter sx dupl sr
if ug eq 0 enter sy dupl sr
set k 1
label elp
if ub eq 0 then
enter w@K vect z A @K H@{C}R A @K C2@C
enter a@K vect r A @K H@{C}R A @K C2@C
enter x@K vect z A @K H@{C}S A @K C2@C
enter b@K vect r A @K H@{C}S A @K C2@C
else
enter w@K vect z A @K H@{C}1 A @K C2@C
enter a@K vect r A @K H@{C}1 A @K C2@C
endif
if c gt 15 goto eskip
if ug eq 0 then
enter y@K vect z A @K H@{C}X A @K C3@C
enter c@K vect r A @K H@{C}X A @K C3@C
enter z@K vect z A @K H@{C}Y A @K C3@C
enter d@K vect r A @K H@{C}Y A @K C3@C
else
enter y@K vect z A @K H@{C}1 A @K C3@C
enter c@K vect r A @K H@{C}1 A @K C3@C
endif
label eskip
incr k by 1
if k le @NL goto elp

! OPEN FILES, FILL SERIES VIA TRAJ

open read unit 1 file name popclipidsnowater.dcd

! LOOP OVER LIPIDS; NORM, CALC SCD, ACCUM
set k 1
label clp
mantim w@K ratio a@K
mantim w@K squa
mantim w@K mult 1.5
mantim w@K shift -0.5
mantim sr add w@K
if ub eq 0 then
mantim x@K ratio b@K
mantim x@K squa
mantim x@K mult 1.5
mantim x@K shift -0.5
mantim ss add x@K
endif
if c gt 15 goto mskip
mantim y@K ratio c@K
mantim y@K squa
mantim y@K mult 1.5
mantim y@K shift -0.5
mantim sx add y@K
if ug eq 0 then
mantim z@K ratio d@K
mantim z@K squa
mantim z@K mult 1.5
mantim z@K shift -0.5
mantim sy add z@K
endif
label mskip
incr k by 1
if k le @NL goto clp

! use final mantime to set AVER, FLUC for output
open write unit 1 card name cd@{C}.txt
echu 1
mantim sr divi @NL
echo sr@C ?AVER ?FLUC
if ub eq 0 then
mantim ss divi @NL
echo ss@C ?AVER ?FLUC
endif
if c gt 15 goto oskip
mantim sx divi @NL
echo sx@C ?AVER ?FLUC
if ug eq 0 then
mantim sy divi @NL
echo sy@C ?AVER ?FLUC
endif
label oskip
close unit 1

I don't see any problems with the input; however, the error message suggests there are non-printing control characters present in the file. While one possibility for that is accidental file corruption, sometimes this can occur for files which have been edited on a Windows or Mac OS X machine and then transferred to a Linux machine to run CHARMM. Note that TAB characters are not allowed in CHARMM input scripts. You may need to consult with a local computer guru to help you with this.

One thing you can try is to make a new file by using a graphic cut-and-paste from the old file to a new file, using a text editor on the Linux system. In principle, this should not copy any non-printing characters which may be causing trouble.

OK, I fixed that error, but I still have results 0.
Here is a sample of the output:

NONBOND OPTION FLAGS:
ELEC VDW ATOMs CDIElec FSHIft VATOm VFSWIt
BYGRoup NOEXtnd NOEWald
CUTNB = 14.000 CTEXNB =999.000 CTONNB = 10.000 CTOFNB = 12.000
CGONNB = 0.000 CGOFNB = 10.000
WMIN = 1.500 WRNMXD = 0.500 E14FAC = 1.000 EPS = 1.000
NBXMOD = 5
There are 0 atom pairs and 0 atom exclusions.
There are 0 group pairs and 0 group exclusions.
with mode 5 found 106586 exclusions and 97544 interactions(1-4)
found 24112 group exclusions.
Generating nonbond list with Exclusion mode = 5
== PRIMARY == SPACE FOR 10550403 ATOM PAIRS AND 0 GROUP PAIRS
== PRIMARY == SPACE FOR 15825624 ATOM PAIRS AND 0 GROUP PAIRS

General atom nonbond list generation found:
15104476 ATOM PAIRS WERE FOUND FOR ATOM LIST
823754 GROUP PAIRS REQUIRED ATOM SEARCHES


CORREL> enter sr zero

CORREL> if ub eq 0 enter ss dupl sr
Comparing "0" and "0".
IF test evaluated as true. Performing command

CORREL> enter sx dupl sr

CORREL> if ug eq 0 enter sy dupl sr
Comparing "0" and "0".
IF test evaluated as true. Performing command

CORREL>

CORREL> set k 1
Parameter: K <- "1"

CORREL> label elp

CORREL> if ub eq 0 then
Comparing "0" and "0".
IF test evaluated as true. Performing command

CORREL> enter w@K vect z A @K H@{C}R A @K C2@C
Parameter: K -> "1"
Parameter: K -> "1"
Parameter: C -> "2"
Parameter: K -> "1"
Parameter: C -> "2"

CORREL> enter a@K vect r A @K H@{C}R A @K C2@C
Parameter: K -> "1"
Parameter: K -> "1"
Parameter: C -> "2"
Parameter: K -> "1"
Parameter: C -> "2"

CORREL> enter x@K vect z A @K H@{C}S A @K C2@C
Parameter: K -> "1"
Parameter: K -> "1"
Parameter: C -> "2"
Parameter: K -> "1"
Parameter: C -> "2"

CORREL> enter b@K vect r A @K H@{C}S A @K C2@C
Parameter: K -> "1"
Parameter: K -> "1"
Parameter: C -> "2"
Parameter: K -> "1"
Parameter: C -> "2"

CORREL> else
Skip commands until next ENDIF

CORREL> if c gt 15 goto eskip





CORREL> ! LOOP OVER LIPIDS; NORM, CALC SCD, ACCUM
CORREL> set k 1
Parameter: K <- "1"

CORREL> label clp

CORREL> mantim w@K ratio a@K ! DIVIDE VECTOR Z BY R (z OF UNIT VECTOR)
Parameter: K -> "1"
Parameter: K -> "1"
NSER: 1
NAMES: W1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim w@K squa ! z**2
Parameter: K -> "1"
NSER: 1
NAMES: W1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim w@K mult 1.5 ! * 3/2
Parameter: K -> "1"
NSER: 1
NAMES: W1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim w@K shift -0.5 ! - 1/2
Parameter: K -> "1"
NSER: 1
NAMES: W1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim sr add w@K
Parameter: K -> "1"
NSER: 1
NAMES: SR
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

CORREL> if ub eq 0 then
Comparing "0" and "0".
IF test evaluated as true. Performing command

CORREL> mantim x@K ratio b@K
Parameter: K -> "1"
Parameter: K -> "1"
NSER: 1
NAMES: X1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim x@K squa
Parameter: K -> "1"
NSER: 1
NAMES: X1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim x@K mult 1.5
Parameter: K -> "1"
NSER: 1
NAMES: X1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim x@K shift -0.5
Parameter: K -> "1"
NSER: 1
NAMES: X1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim ss add x@K
Parameter: K -> "1"
NSER: 1
NAMES: SS
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

CORREL> endif

CORREL> if c gt 15 goto mskip
Comparing "2" and "15".
IF test evaluated as false. Skipping command

CORREL> mantim y@K ratio c@K
Parameter: K -> "1"
Parameter: K -> "1"
NSER: 1
NAMES: Y1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim y@K squa
Parameter: K -> "1"
NSER: 1
NAMES: Y1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim y@K mult 1.5
Parameter: K -> "1"
NSER: 1
NAMES: Y1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim y@K shift -0.5
Parameter: K -> "1"
NSER: 1
NAMES: Y1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim sx add y@K
Parameter: K -> "1"
NSER: 1
NAMES: SX
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

CORREL> if ug eq 0 then
Comparing "0" and "0".
IF test evaluated as true. Performing command

CORREL> mantim z@K ratio d@K
Parameter: K -> "1"
Parameter: K -> "1"
NSER: 1
NAMES: Z1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim z@K squa
Parameter: K -> "1"
NSER: 1
NAMES: Z1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim z@K mult 1.5
Parameter: K -> "1"
NSER: 1
NAMES: Z1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim z@K shift -0.5
Parameter: K -> "1"
NSER: 1
NAMES: Z1
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: VECT
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 3
VALUE: 0.000000

CORREL> mantim sy add z@K
Parameter: K -> "1"
NSER: 1
NAMES: SY
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

CORREL> endif

CORREL> label mskip

CORREL> incr k by 1
Parameter: K <- "2"

CORREL> if k le @NL goto clp
Parameter: NL -> "274"
Comparing "2" and "274".
IF test evaluated as true. Performing command



CORREL> ! use final mantime to set AVER, FLUC for output
CORREL> open write unit 1 card name cd@{C}.txt
Parameter: C -> "2"
OPNLGU> Unit already open. The old file will be closed first.
VCLOSE: Closing unit 1 with status "KEEP"
VOPEN> Attempting to open::cd2.txt::
OPNLGU> Unit 1 opened for WRITE access to cd2.txt

CORREL> echu 1

CORREL> mantim sr divi @NL
Parameter: NL -> "274"
NSER: 1
NAMES: SR
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

CORREL> echo sr@C ?AVER ?FLUC
Parameter: C -> "2"
RDCMND substituted energy or value "?AVER" to "0"
RDCMND substituted energy or value "?FLUC" to "0"

CORREL> if ub eq 0 then
Comparing "0" and "0".
IF test evaluated as true. Performing command

CORREL> mantim ss divi @NL
Parameter: NL -> "274"
NSER: 1
NAMES: SS
TOTALS: 0
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

CORREL> echo ss@C ?AVER ?FLUC
Parameter: C -> "2"
RDCMND substituted energy or value "?AVER" to "0"
RDCMND substituted energy or value "?FLUC" to "0"

CORREL> endif

CORREL> if c gt 15 goto oskip

Why do I keep getting 0 results?

SR2 0 0
SS2 0 0
SX2 0 0
SY2 0 0

Please review the READ BEFORE POSTING guidelines for suggestions on the best way to post output logs.

It looks like you may have omitted the all-important TRAJ command, which actually reads the trajectory file(s) and fills the requested time series with data. The comments in the sample inputs are usually important, such as this one:

! OPEN FILES, FILL SERIES VIA TRAJ

I didn't ommit. In the script I listed above I included TRAJ command as:

! OPEN FILES, FILL SERIES VIA TRAJ
open read unit 1 file name popclipidsnowater.dcd

I have one file only, popclipidsnowater.dcd. What else should I include?

Also, in the output, I got no errors in opening the trajectory file:
CORREL> ! OPEN FILES, FILL SERIES VIA TRAJ
CORREL> open read unit 1 file name popclipidsnowater.dcd
OPNLGU> Unit already open. The old file will be closed first.
VCLOSE: Closing unit 1 with status "KEEP"
VOPEN> Attempting to open::popclipidsnowater.dcd::
OPNLGU> Unit 1 opened for READONLY access to popclipidsnowater.dcd

CORREL>


Thank you,

The above only opens the file; look more closely at the original script.

OK, this is the original script:
! OPEN FILES, FILL SERIES VIA TRAJ
set k = @FF
set u = 8
label opnlp
open read unit @U file name @S/dyn@K.trj
incr u by 1
incr k by 1
if k le @LF goto opnlp
traj firstu 8 nunit @N

I see here a loop to merge several trajectory files. In my case I don't have several trajectory file to merge. I have only one file, popclipidsnowater.dcd.
I don't see other command besides "open read unit @U file name @S/dyn@K.trj"

What am I missing?

Thank you,

There is no merge of several trajectory files. The loop simply opens the files ("connects the name to a unit number"). The actual reading of the file(s) is done withe the "traj firstu ..." command. What did you think this command was doing?

Since you're using a developer's version of CHARMM, there has to be someone in your institution who can help you answer all these basic questions.

What causes value 0 in fluctuation of order parameter?

I get all fluctuations as 0, starting like this:

TITLE> *
***** WARNING ***** BEGIN= 0 Was not specified. It has been set to: 0

1 CORD RECORDS READ FROM 1 UNITS STARTING WITH UNIT 1
RUNNING FROM STEP 0 TO 0 SKIPPING 1 STEPS BETWEEN RECORDS
Time step was 1.000000 AKMA time units.
1 Series "SR " Average = 0.000000 rms Fluctuation = 0.000000
2 Series "SS " Average = 0.000000 rms Fluctuation = 0.000000
3 Series "SX " Average = 0.000000 rms Fluctuation = 0.000000
4 Series "SY " Average = 0.000000 rms Fluctuation = 0.000000
5 Series "W1 " Average = -0.524498 rms Fluctuation = 0.000000
6 Series "A1 " Average = 1.110998 rms Fluctuation = 0.000000
7 Series "X1 " Average = -0.133732 rms Fluctuation = 0.000000
8 Series "B1 " Average = 1.111001 rms Fluctuation = 0.000000
9 Series "Y1 " Average = 0.705927 rms Fluctuation = 0.000000
10 Series "C1 " Average = 1.111002 rms Fluctuation = 0.000000
11 Series "Z1 " Average = -0.664016 rms Fluctuation = 0.000000

Why is that?

One more thing, just to be sure: I suppose I have to calculate an average between SR and SS and between SX and SY? (That in order to plot order parameter of each tail vs. carbon number).
Thank you,

The fluctuations are 0.0 because you've only read one frame. The script was written to analyze an ensemble of thousands of frames.

The averages are computed in the script and written to the output file.

These really are very, very basic questions that you should be discussing with others in your lab rather than posting here.

In which way I specified reading one frame only and how can I correct that?

Thank you,

What if you look across departmental boundaries?

correl.doc provides the syntax for its traj subcommand; see also dynamc.doc.

We do not know how you specified only one frame, because you have not shown us the traj command you used. Does your trajectory contain more than one frame?

I should note that the script computes averages for a given C-H vector over all lipids and time points, but does not average the two values for the two vectors on the same C atom; you would need to that. They are not averaged because the equivalence of the two values is an indication of convergence; for a well converged analysis interval, the values should agree fairly closely.

Hi, can we use above code to calculate the order parameters of PIP2? It looks like it also has the specified carbon atoms.