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

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?