CHARMM Development Project Forums CHARMM Community Script Archive merge-orient-trajectory.inp

(600*(32/2))/250

If I plug into that formula frame 1468 (the last frame of the merge traj, which should be the last trajectory), (1468*(32/2))/250=94 I should expect 110=150-40=(final traj in merge)-(initial traj in merge)

Not 1438 frames, as indicated above?

No. of files = 1 + (last-first), so in principle you've merged 111 files. The expected total from the merge, assuming no errors should be

111*250*(2/32) or 1734 frames

It doesn't add up; you're going to have to actually look at the output log files.

My incorrect "nunit" is the cause. Thanks.

My reference structure has one atom less than the trajectory. The reference and the trajectory structures are different protein where it differs with one extra ion in it.

I need to merge 150 .crd files (883 atoms) into a single .dcd file after least square fitting it with a reference structure (882 atoms).

i was using the modified merge-orient-trajectory in the basic script. The error is number of atoms mismatch. I need to write all the atoms present in the .crd to the dcd. So , please help me.

Is it possible

MY code is below

bomlev -1
set dcd = ./dcd
set fid = cg.polar
set refdir = /home/nixon/phd/nw_flucts/proa/

read rtf card name top.inp
read para card name par.inp
read psf card name @refdir/@fid.psf
read coor card name @refdir/@fid.crd

! put reference coordinates into comparison
coor copy comp

! New, merged trajectory
open write unit 160 file name merged.dcd
traj iwrite 160 nwrite 1 nfile 149 skip 1
* trajectory of merged crd files
*

set 1 1
set 2 11
LABEL LOOP1
open read unit @2 form name @dcd/cg.@1.crd
read coor card unit @2
close unit @2
traj write
incr 1 BY 1
incr 2 by 1
IF 1 LT 150 GOTO LOOP1

define mysel sele segid PROA .and. .not. resname A299 end

merge first 11 nunit 149 output 160 skip 1 - ! one frame/ps is all we need
sele segid PROB end orient mass sele mysel end

stop

[quote=Nixon]MY code is below

bomlev -1
set dcd = ./dcd
set fid = cg.polar
set refdir = /home/nixon/phd/nw_flucts/proa/

read rtf card name top.inp
read para card name par.inp
read psf card name @refdir/@fid.psf
read coor card name @refdir/@fid.crd

! put reference coordinates into comparison
coor copy comp

! New, merged trajectory
open write unit 160 file name merged.dcd
traj iwrite 160 nwrite 1 nfile 149 skip 1
* trajectory of merged crd files
*

set 1 1
set 2 11
LABEL LOOP1
open read unit @2 form name @dcd/cg.@1.crd
read coor card unit @2
close unit @2
traj write
incr 1 BY 1
incr 2 by 1
IF 1 LT 150 GOTO LOOP1

define mysel sele segid PROA .and. .not. resname A299 end

merge first 11 nunit 149 output 160 skip 1 - ! one frame/ps is all we need
sele segid PROB end orient mass sele mysel end

stop

I get this error

At line 2111 of file /c37b2/source/dynamc/cvio.src (unit = 11, file = 'fort.11')
Fortran runtime error: End of file

And this is the final part of my log file output

CHARMM> IF 1 LT 150 GOTO LOOP1
Comparing "149" and "150".
IF test evaluated as true. Performing command

CHARMM> open read unit @2 form name @dcd/cg.@1.crd
Parameter: 2 -> "159"
Parameter: DCD -> "./DCD"
Parameter: 1 -> "149"
VOPEN> Attempting to open::./dcd/cg.149.crd::
OPNLGU> Unit159 opened for READONLY access to ./dcd/cg.149.crd

CHARMM> read coor card unit @2
Parameter: 2 -> "159"
SPATIAL COORDINATES BEING READ FROM UNIT159
TITLE> * FRAME 4000 FROM /HOME/NIXON/PHD/NW_FLUCTS/3CM5/MONO_3CM5/POLAR/DCD/CG.149.DCD
TITLE> *

** WARNING ** For atom in coordinate file, the residue number is out of range, and is thus ignored:

IRES= 300 RESNAME= A300 TYPE= A883
*** LEVEL 1 WARNING *** BOMLEV IS -1

** WARNING ** Coordinates were overwritten for 882 atoms.

** MESSAGE ** 1 atoms in coordinate file were outside the specified sequence range.
*** LEVEL 2 WARNING *** BOMLEV IS -1

CHARMM> close unit @2
Parameter: 2 -> "159"
VCLOSE: Closing unit 159 with status "KEEP"

CHARMM> traj write

149 COORDINATE SETS STARTING FROM
STEP NO 1 FOR EVERY 1 STEPS
WRITTEN ON UNIT 160

VCLOSE: Closing unit 160 with status "KEEP"

CHARMM> incr 1 BY 1
Parameter: 1 <- "150"

CHARMM> incr 2 by 1
Parameter: 2 <- "160"

CHARMM> IF 1 LT 150 GOTO LOOP1
Comparing "150" and "150".
IF test evaluated as false. Skipping command

CHARMM>

CHARMM> define mysel sele segid 3cg7 .and. .not. resname A299 end
SELRPN> 881 atoms have been selected out of 882

CHARMM>

CHARMM> merge first 11 nunit 149 output 160 skip 1 - ! one frame/ps is all we need
CHARMM> sele segid 3cm5 end orient mass sele mysel end

merge first 11
unit 11 is not a trajetctory - as far as I can see this is one of the 150 coordinate files.

It looks like you may be intermixing TRAJ I/O commands with MERGE commands.

Note also that SKIP values are in integration steps, not time units.

When reading trajectory files, the number and ordering of atoms must match exactly what is in the PSF. The MERGE command does allow the writing of a subset of atoms, but requires another trajectory as input data.