CHARMM Development Project
This is untested, but my concern is the enter dist command.
I am selecting multiple atoms hoping for it to recognize that I want the center of mass of the selection. THanks

Code:
stream toppar.str

read psf card name "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.psf"
OPEN UNIT 1 READ CARD NAME "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.crd"
READ COOR UNIT 1 CARD
CLOSE UNIT 1

coor copy comp

OPEN READ UNIT 100 FILE NAME "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.dcd"

Set maxmol 54

correl maxtime 1000000 maxseries 100

Label mol
!                   segid resid    type segid resid     type
enter v@molnum dist DPC   @molnum  P    DPC   1:@maxmol *
If molnum LT @maxmol GOTO mol

TRAJECTORY FIRSTU 100 NUNIT 1

!add add timeseries

Set molnum 2

Label mol
mantime v1 ADD v@molnum
If molnum LT @maxmol GOTO mol

mantime v1 divi maxmol

I don't think that will work-- I haven't checked the Fortran source code, but I believe that only atom selection supports wildcards (*) and specifying ranges (:). Also, correl.doc implies that the DIST time series expects only two atoms.

You probably need to use ATOM time series to get the P atom coordinates and the center of mass (COM) coordinates, then use MANTIME commands to subtract the COM from the P atom coordinates, and then compute the vector length.

In order to compute the vector length within CORREL, it may be easiest (although cumbersome) to specify the individual x,y,z components instead of the ordered xyz triple, e.g.

enter comx atom x sele segid DPC end mass
enter comy atom y sele segid DPC end mass
enter comz atom z sele segid DPC end mass


Here, the atom selection is advantageous, as the ATOM series will compute the average coordinate(s) over the list of atoms provided. (Not all time series types will do this.)
Although there are no errors, the resulting time series is smaller than expected; all values are less than 1.
I must have done something wrong in recreating the distance formula.



Code:
stream toppar.str

read psf card name "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.psf"
OPEN UNIT 1 READ CARD NAME "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.crd"
READ COOR UNIT 1 CARD
CLOSE UNIT 1

coor copy comp

open unit 11 write form name distpcom.dat
write title unit 11
* time distpcom
*

OPEN READ UNIT 100 FILE NAME "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.dcd"

correl maxtime 1000000 maxseries 100

enter comx atom x sele segid DPC end mass
enter comy atom y sele segid DPC end mass
enter comz atom z sele segid DPC end mass
enter px atom x sele type P end mass
enter py atom y sele type P end mass
enter pz atom z sele type P end mass

TRAJECTORY FIRSTU 100 NUNIT 1 stop 74002000


!make negative
mantime comx mult -1
mantime comy mult -1
mantime comz mult -1

!subtract by adding the negative
mantime px add comx
mantime py add comy
mantime pz add comz

!square the difference
mantime px ipower 2
mantime py ipower 2
mantime pz ipower 2

!sum
mantime px add py
mantime px add pz

!square root
mantime px sqrt

write px dumb time unit 11
*hi
*
end

The px,py,pz time series have the average position of all P atoms, which is probably close to the COM for an intact micelle. You need to get the vector length to each individual P atom.
I made no changes to the loading of the parameter/forcefield/topology files, but I got an error in that portion:

Code:
[login-2:12549] *** Process received signal ***
[login-2:12549] Signal: Segmentation fault (11)
[login-2:12549] Signal code: Address not mapped (1)
[login-2:12549] Failing at address: 0xbadda73c38
[login-2:12549] [ 0] /lib64/libpthread.so.0 [0x2b5fb68a0ca0]
[login-2:12549] [ 1] /cell_root/software/openmpi/1.2.5/gnu/sys/lib/libopen-pal.so.0(_int_malloc+0xb39) [0x2b5fb59c14b9]
[login-2:12549] [ 2] /cell_root/software/openmpi/1.2.5/gnu/sys/lib/libopen-pal.so.0(malloc+0x88) [0x2b5fb59c2a38]
[login-2:12549] [ 3] /homes/jbklauda/charmm/c35b5/exec/gnu/g77-xlg(veheap_+0xd8) [0x6b0838]
[login-2:12549] [ 4] /homes/jbklauda/charmm/c35b5/exec/gnu/g77-xlg(allhp_+0x28d) [0xb7109d]
[login-2:12549] [ 5] /homes/jbklauda/charmm/c35b5/exec/gnu/g77-xlg(correl_+0x43e) [0x45cefe]
[login-2:12549] [ 6] /homes/jbklauda/charmm/c35b5/exec/gnu/g77-xlg(maincom_+0x319) [0x437789]
[login-2:12549] [ 7] /homes/jbklauda/charmm/c35b5/exec/gnu/g77-xlg(MAIN__+0xad5) [0x43aa55]
[login-2:12549] [ 8] /homes/jbklauda/charmm/c35b5/exec/gnu/g77-xlg(main+0xe) [0xc1f09e]
[login-2:12549] [ 9] /lib64/libc.so.6(__libc_start_main+0xf4) [0x2b5fb6acb9c4]
[login-2:12549] [10] /homes/jbklauda/charmm/c35b5/exec/gnu/g77-xlg [0x4348c9]
[login-2:12549] *** End of error message ***
Segmentation fault


My updated script that I tried to run:

Code:
stream toppar.str

read psf card name "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.psf"
OPEN UNIT 1 READ CARD NAME "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.crd"
READ COOR UNIT 1 CARD
CLOSE UNIT 1

coor copy comp

open unit 11 write form name distpcom.dat
write title unit 11
* time distpcom
*

OPEN READ UNIT 100 FILE NAME "/export/lustre_1/mallsopp/dpc/radius/setup/micelle-tmp-sim1.dcd"

correl maxtime 1000000 maxseries 100000

enter comx atom x sele segid DPC end mass
enter comy atom y sele segid DPC end mass
enter comz atom z sele segid DPC end mass

set num 1
label loop

enter px@num atom x sele resid @num .and. type P end mass
enter py@num atom y sele resid @num .and. type P end mass
enter pz@num atom z sele resid @num .and. type P end mass

incr num by 1
if num lt @nres goto loop

TRAJECTORY FIRSTU 100 NUNIT 1 stop 74002000


!make negative
mantime comx mult -1
mantime comy mult -1
mantime comz mult -1

set num 1
label loop2

!subtract by adding the negative
mantime px@num add comx
mantime py@num add comy
mantime pz@num add comz

!square the difference
mantime px@num ipower 2
mantime py@num ipower 2
mantime pz@num ipower 2

!sum
mantime px@num add py@num
mantime px@num add pz@num

!square root
mantime px@num sqrt

incr num by 1
if @num lt @nres goto loop2

!now average all px@num 's

set num 2
label loop3
mantime px1 add px@num
incr num by 1
if num lt @nres goto loop3

mantime px1 divi @nres

write px1 dumb time unit 11
*hi
*
end
It can be difficult to diagnose the cause of seg fault errors, esp. out of context like this. However, the call trace indicates the failure may be related to a problem with HEAP expansion, an issue with legacy versions such as c35b5, released prior to the Fortran95 conversion of CHARMM. (I tend to use a non-parallel executable for analysis; the call trace includes openmpi libs).

Note that your MAXTIME and MAXSER values are requesting 745 GB of memory, a likely cause of a HEAP expansion seg fault. Their product determines the memory allocation, which I computed as (8*MAXTIME*MAXSER)/(1024*1024*1024).

It's also possible that MAXATOM may be too small, but the MAXTIME and MAXSER values need to fixed first. The value for MAXSER can be computed via

calc mxs = 3 * ( @NRES + 1 )

Note that

sele resid @num .and. type P end mass

can be done more simply as an atom-spec:

DPC @NUM P

For a single atom, there's no need for the MASS keyword.
Thanks very much. I got some reasonable results.

Unexpectedly, I had to set nres.
Subst.doc led me to believe
it would automatically be set after the psf was read.
it is ?NRES, not @NRES (? means a value set internally, and @ means a value set by the user). subst.doc deals with the first kind.
It seems that except for the COM, all time series are zero.
Check your output log for messages.
I am getting a good result.

Since the psf contains water and DPC, then I should not use ?nres because it include waters.

However the timeseries is still influenced by the intermittent breakup by reaching the boundary of the unit cell.
But it is easier to exclude those values.
The problem is that the COM coordinates will also be a bit off because of lipids hopping across the unit cell due to image centering.
Yes, but if I can distinguish the snapshots that this is occurring and not include them it will be alright, right?
I can't upload the excel format, but heres a screenshot.

[img:right]http://imgur.com/D8DSrME[/img]
If you only include coordinate sets where the micelle is intact, then I believe you should get a reasonable estimate of the average radial distance of each P atom from the center of mass.
© CHARMM forums