|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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
Rick Venable computational chemist
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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
|
|
|
|
Joined: Jul 2004
Posts: 165
Forum Member
|
Forum Member
Joined: Jul 2004
Posts: 165 |
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.
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Jul 2004
Posts: 113
Forum Member
|
Forum Member
Joined: Jul 2004
Posts: 113 |
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.
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Jul 2004
Posts: 113
Forum Member
|
Forum Member
Joined: Jul 2004
Posts: 113 |
Which one is the most important factor among the above three ? Is it necessary to check the convergence of the order parameter ?
Thanks.
|
|
|
|
Joined: Jan 2005
Posts: 86
Forum Member
|
Forum Member
Joined: Jan 2005
Posts: 86 |
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.
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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
Last edited by rmv; 07/28/14 10:34 PM. Reason: no. of files update
Rick Venable computational chemist
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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
Rick Venable computational chemist
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
Hello,
What does "Error in parameter substitution" mean and how can it be fixed?
Thank you, Mihaela
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
The value for @C is passed as a command line argument; see the very first post in this thread for an example.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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?
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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?
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
There are probably messages in the output file which give some indication of the problem.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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?
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
I usually keep the water.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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?
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
Find the clues in your output file.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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?
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
Why do I keep getting 0 results?
SR2 0 0 SS2 0 0 SX2 0 0 SY2 0 0
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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,
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
The above only opens the file; look more closely at the original script.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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,
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
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?
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Dec 2005
Posts: 1,535
Forum Member
|
Forum Member
Joined: Dec 2005
Posts: 1,535 |
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.
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
What causes value 0 in fluctuation of order parameter?
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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?
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
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,
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Feb 2013
Posts: 25
Forum Member
|
Forum Member
Joined: Feb 2013
Posts: 25 |
In which way I specified reading one frame only and how can I correct that?
Thank you,
Last edited by Mih; 10/09/14 09:12 PM.
|
|
|
|
Joined: Dec 2005
Posts: 1,535
Forum Member
|
Forum Member
Joined: Dec 2005
Posts: 1,535 |
|
|
|
|
Joined: Sep 2003
Posts: 4,881 Likes: 12
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 4,881 Likes: 12 |
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?
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Sep 2003
Posts: 8,650 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,650 Likes: 26 |
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.
Rick Venable computational chemist
|
|
|
|
Joined: Aug 2009
Posts: 139
Forum Member
|
Forum Member
Joined: Aug 2009
Posts: 139 |
Hi, can we use above code to calculate the order parameters of PIP2? It looks like it also has the specified carbon atoms.
|
|
|
|
|