CHARMM Development Project
lipid C-D order parameters using CORREL - 05/26/05 01:36 AM
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
Re: lipid C-D order parameters using CORREL - 05/26/05 02:01 AM

`* axial avg order parameter for ethanolamine C:H vectors; @N files*bomblev -1! READ STD FILES FOR PROJECT; RTF, PARAM, PSF, COORstream rtfprm.strstream psfcrd.str! COMPUTE CORREL LIMITScalc mxa 80 * 6 * 4calc mxs 10 + 80*2 * 4calc mxt 200 * @N! INVOKE CORRELcorrel maxtim @MXT maxa @MXA maxs @MXS! ACCUMULATE AVERAGES IN THESE TIME SERIES; DECLARE IN THIS ORDERenter sr zeroenter ss dupl srenter sx dupl srenter sy dupl sr! LOOP OVER THE LIPIDS; USE ONLY Z AND R OF VECTORset k 1label elpenter w@K vect z L @K H11A  L @K C11enter a@K vect r L @K H11A  L @K C11enter x@K vect z L @K H11B  L @K C11enter b@K vect r L @K H11B  L @K C11enter y@K vect z L @K H12A  L @K C12enter c@K vect r L @K H12A  L @K C12enter z@K vect z L @K H12B  L @K C12enter d@K vect r L @K H12B  L @K C12incr k by 1if k le 80 goto elp! START FROM UNIT 8; ALLOWS UP TO CA. 90 FILESset k 1label klpcalc u @K + 7open unit @U file read name dyn@K.trjincr k by 1if k .le. @N goto klptraj firstu 8 nunit @N! LOOP OVER 80 LIPIDS; NORM, CALC SCD, ACCUMcalc r80 1. / 80.set k 1label clpmantim w@K ratio a@K   ! DIVIDE VECTOR Z BY R (z OF UNIT VECTOR)mantim w@K squa        ! z**2mantim w@K mult 1.5    ! * 3/2mantim w@K shift -0.5  ! - 1/2mantim w@K mult @R80   ! * 1/Nlipidmantim sr add w@K      ! ACCUMLATE, H11Amantim x@K ratio b@Kmantim x@K squamantim x@K mult 1.5mantim x@K shift -0.5mantim x@K mult @R80mantim ss add x@K      ! ACCUMLATE, H11Bmantim y@K ratio c@Kmantim y@K squamantim y@K mult 1.5mantim y@K shift -0.5mantim y@K mult @R80mantim sx add y@K      ! ACCUMLATE, H12Amantim z@K ratio d@Kmantim z@K squamantim z@K mult 1.5mantim z@K shift -0.5mantim z@K mult @R80mantim sy add z@K      ! ACCUMLATE, H12Bincr k by 1if k le 80 goto clp! LOOP OVER LIPIDS ENDS; COMBINE THE FOUR TIME SERIESedit sr veccod 4 delta 0.001 skip 1000 offset 1.! STASH DATA IN A SUBDIR NAMED scd; 5 COLS, TIME IN 1st COLopen write unit 1 card name scd/eth.datwrite sr unit 1 dumb time* dumb*endstop`
Re: lipid C-D order parameters using CORREL - 05/26/05 02:06 AM
`* the axial avg order parameter for glycerol C:H vectors; @N files*bomblev -2stream rtfprm.strstream psfcrd.strcalc mxa 80 * 6 * 5calc mxs 10 + 80*2 * 5calc mxt 200 * @Ncorrel maxtim @MXT maxa @MXA maxs @MXSenter sa zeroenter sb dupl saenter ss dupl saenter sx dupl saenter sy dupl saset k 1label elpenter v@K vect z L @K HS  L @K C2enter r@K vect r L @K HS  L @K C2enter w@K vect z L @K HA  L @K C1enter a@K vect r L @K HA  L @K C1enter x@K vect z L @K HB  L @K C1enter b@K vect r L @K HB  L @K C1enter y@K vect z L @K HX  L @K C3enter c@K vect r L @K HX  L @K C3enter z@K vect z L @K HY  L @K C3enter d@K vect r L @K HY  L @K C3incr k by 1if k le 80 goto elpset k 1label klpcalc u @K + 7open unit @U file read name dyn@K.trjincr k by 1if k .le. @N goto klptraj firstu 8 nunit @N! LOOP OVER 80 LIPIDS; NORM, CALC SCD, ACCUMcalc r80 1. / 80.set k 1label clpmantim w@K ratio a@Kmantim w@K squamantim w@K mult 1.5mantim w@K shift -0.5mantim w@K mult @R80mantim sa add w@Kmantim x@K ratio b@Kmantim x@K squamantim x@K mult 1.5mantim x@K shift -0.5mantim x@K mult @R80mantim sb add x@Kmantim y@K ratio c@Kmantim y@K squamantim y@K mult 1.5mantim y@K shift -0.5mantim y@K mult @R80mantim sx add y@Kmantim z@K ratio d@Kmantim z@K squamantim z@K mult 1.5mantim z@K shift -0.5mantim z@K mult @R80mantim sy add z@Kmantim v@K ratio r@Kmantim v@K squamantim v@K mult 1.5mantim v@K shift -0.5mantim v@K mult @R80mantim ss add v@Kincr k by 1if k le 80 goto clpedit sa veccod 5 delta 0.001 skip 1000 offset 1.open write unit 1 card name scd/gly.datwrite sa unit 1 dumb time* dumb*end`
Re: lipid C-D order parameters using CORREL - 05/26/05 02:14 AM
`* the axial avg order parameter for an input C no.; @N files* additional required arg C: chain carbon number @C*bomblev -2stream rtfprm.strstream psfcrd.strcalc mxa 80 * 6 * 4calc mxs 10 + 80*2 * 4calc mxt 200 * @Ncorrel maxtim @MXT maxa @MXA maxs @MXSenter sr zero enter ss dupl srenter sx dupl srenter sy dupl srset k 1! TRAILING LETTERS IN THE 'H' ATOM NAMES REQUIRE @{C}, NOT @Clabel elpenter w@K vect z L @K H@{C}R  L @K C2@Center a@K vect r L @K H@{C}R  L @K C2@Center x@K vect z L @K H@{C}S  L @K C2@Center b@K vect r L @K H@{C}S  L @K C2@Center y@K vect z L @K H@{C}X  L @K C3@Center c@K vect r L @K H@{C}X  L @K C3@Center z@K vect z L @K H@{C}Y  L @K C3@Center d@K vect r L @K H@{C}Y  L @K C3@Cincr k by 1if k le 80 goto elpset k 1label klpcalc u @K + 7open unit @U file read name dyn@K.trjincr k by 1if k .le. @N goto klptraj firstu 8 nunit @N! LOOP OVER LIPIDS; NORM, CALC SCD, ACCUMcalc r80 1. / 80.set k 1label clpmantim w@K ratio a@Kmantim w@K squamantim w@K mult 1.5mantim w@K shift -0.5mantim w@K mult @R80mantim sr add w@Kmantim x@K ratio b@Kmantim x@K squamantim x@K mult 1.5mantim x@K shift -0.5mantim x@K mult @R80mantim ss add x@Kmantim y@K ratio c@Kmantim y@K squamantim y@K mult 1.5mantim y@K shift -0.5mantim y@K mult @R80mantim sx add y@Kmantim z@K ratio d@Kmantim z@K squamantim z@K mult 1.5mantim z@K shift -0.5mantim z@K mult @R80mantim sy add z@Kincr k by 1if k le 80 goto clpedit 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.datwrite sr unit 1 dumb time* dumb*endstop`
Re: lipid C-D order parameters using CORREL - 08/12/05 05:40 AM
Dear Dr Venable,

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.
Re: lipid C-D order parameters using CORREL - 08/12/05 03:07 PM
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.
Re: lipid C-D order parameters using CORREL - 07/31/06 03:42 PM
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.
Re: lipid C-D order parameters using CORREL - 07/31/06 06:10 PM
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.
Re: lipid C-D order parameters using CORREL - 07/31/06 07:41 PM
Which one is the most important factor among the above three ?
Is it necessary to check the convergence of the order parameter ?

Thanks.
Re: lipid C-D order parameters using CORREL - 08/01/06 06:38 AM
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.

regards,
MzZt.
Re: lipid C-D order parameters using CORREL - 08/01/06 03:19 PM
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 FILESstream rtfprm.str! READ PSF, COORstream psfcrd.str! AND SETUP CRYSTAL (IMPORTANT!!!)stream cryst.str! COMPUTE LAST OF SETcalc ilst @J + 9! START FROM UNIT 21set iu 21! INPUT FILE OPEN LOOPlabel flpopen unit @IU read file name dyn@J.trjincr iu by 1incr j by 1if j .le. @ILST goto flp! OPEN THE OUTPUT FILE AS UNIT 19open unit 19 write file name /scratch/dyn@K.trj! MERGE THE FILES TO CREATE A NEW ONEmerge coor firstu 21 nunit 10 output 19 nfile 1000stop`
Re: lipid C-D order parameters using CORREL - 07/28/14 10:43 PM
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
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
endif
mantim y@K ratio c@K
mantim y@K squa
mantim y@K mult 1.5
mantim y@K shift -0.5
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
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
Re: lipid C-D order parameters using CORREL - 07/28/14 10:53 PM
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
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
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
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
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
Re: lipid C-D order parameters using CORREL - 09/30/14 02:42 PM
Hello,

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

Thank you,
Mihaela
Re: lipid C-D order parameters using CORREL - 09/30/14 03:51 PM
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.
Re: lipid C-D order parameters using CORREL - 09/30/14 07:11 PM
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

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

Thank you,
Mihaela
Re: lipid C-D order parameters using CORREL - 09/30/14 09:00 PM
The value for @C is passed as a command line argument; see the very first post in this thread for an example.
Re: lipid C-D order parameters using CORREL - 10/02/14 03:50 PM
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
Re: lipid C-D order parameters using CORREL - 10/02/14 04:33 PM
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.

Re: lipid C-D order parameters using CORREL - 10/02/14 05:45 PM
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?
Re: lipid C-D order parameters using CORREL - 10/02/14 06:11 PM
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.
Re: lipid C-D order parameters using CORREL - 10/02/14 06:20 PM
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
Re: lipid C-D order parameters using CORREL - 10/02/14 06:21 PM
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
Re: lipid C-D order parameters using CORREL - 10/02/14 07:48 PM
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.
Re: lipid C-D order parameters using CORREL - 10/02/14 09:40 PM
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
Re: lipid C-D order parameters using CORREL - 10/02/14 09:47 PM
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.
Re: lipid C-D order parameters using CORREL - 10/02/14 10:04 PM
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?
Re: lipid C-D order parameters using CORREL - 10/02/14 10:19 PM
There are probably messages in the output file which give some indication of the problem.
Re: lipid C-D order parameters using CORREL - 10/02/14 11:19 PM
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?
Re: lipid C-D order parameters using CORREL - 10/03/14 12:57 AM
I usually keep the water.
Re: lipid C-D order parameters using CORREL - 10/03/14 05:11 PM
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?
Re: lipid C-D order parameters using CORREL - 10/03/14 07:57 PM
Find the clues in your output file.
Re: lipid C-D order parameters using CORREL - 10/03/14 08:41 PM
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::

I think this is the clue, but how do I interpret it?
Re: lipid C-D order parameters using CORREL - 10/04/14 01:22 AM
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.
Re: lipid C-D order parameters using CORREL - 10/04/14 02:45 PM
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

open unit 1 read form name top_all36_lipid.rtf
close unit 1
open unit 1 read form name par_all36_lipid.prm
close unit 1

! Read parameter file in command line

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

! 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
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
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
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
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
Re: lipid C-D order parameters using CORREL - 10/04/14 06:19 PM
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.
Re: lipid C-D order parameters using CORREL - 10/04/14 09:42 PM
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

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

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

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

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::

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
Re: lipid C-D order parameters using CORREL - 10/04/14 09:42 PM
Why do I keep getting 0 results?

SR2 0 0
SS2 0 0
SX2 0 0
SY2 0 0
Re: lipid C-D order parameters using CORREL - 10/05/14 01:49 AM
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
Re: lipid C-D order parameters using CORREL - 10/05/14 06:12 PM
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::

CORREL>

Thank you,
Re: lipid C-D order parameters using CORREL - 10/05/14 09:56 PM
The above only opens the file; look more closely at the original script.
Re: lipid C-D order parameters using CORREL - 10/06/14 12:43 PM
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,
Re: lipid C-D order parameters using CORREL - 10/06/14 03:50 PM
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?
Re: lipid C-D order parameters using CORREL - 10/06/14 06:20 PM
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.
Re: lipid C-D order parameters using CORREL - 10/06/14 11:07 PM
What causes value 0 in fluctuation of order parameter?
Re: lipid C-D order parameters using CORREL - 10/06/14 11:45 PM
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?
Re: lipid C-D order parameters using CORREL - 10/07/14 03:51 AM
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,
Re: lipid C-D order parameters using CORREL - 10/07/14 04:10 AM
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.
Re: lipid C-D order parameters using CORREL - 10/07/14 11:40 AM

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

Thank you,
Re: lipid C-D order parameters using CORREL - 10/07/14 03:27 PM
What if you look across departmental boundaries?
Re: lipid C-D order parameters using CORREL - 10/07/14 03:49 PM