Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
about reading frame from dcd one by one
#27444 05/09/11 10:09 PM
Joined: Mar 2010
Posts: 22
Z
zzlai Offline OP
Forum Member
OP Offline
Forum Member
Z
Joined: Mar 2010
Posts: 22
hi, I have several dcd trajectory files to analyse by my own codes, yet I meet a problem.
my own scripts is:(bash programming)

for ((i=1;i<=$end;i++))
do
charmm <getpdb.inp>getpdb.out (this will produce $i.pdb)
my codes $i.pdb
done


In the above scripts, I would like to produce ONE frame for each loop, so my getpdb.inp is:

!Read topology and parameter
open read card unit 10 name toph19.inp
read rtf card unit 10

open read card unit 20 name param19.inp
read para card unit 20


!Read PSF and initial/reference coordinates
read psf card name protein.psf
read coor pdb resid name protein.pdb
coor copy comp

!Open file unit for trajectory
open read unit 13 file name run1.dcd
open read unit 14 file name run2.dcd
open read unit 15 file name run3.dcd
open read unit 16 file name run4.dcd
open read unit 17 file name run5.dcd

TRAJECTORY FIRSTU 13 NUNIT 5 IREAD BEGIN 500 STOP 500 SKIP 1
TRAJ READ
OPEN WRITE CARD UNIT 12 NAME ${i}.pdb
WRITE COOR PDB UNIT 12
CLOSE UNIT 12

STOP

In each trajectory, there are

READING TRAJECTORY FROM UNIT 13
NUMBER OF COORDINATE SETS IN FILE: 100000
NUMBER OF PREVIOUS DYNAMICS STEPS: 500
FREQUENCY FOR SAVING COORDINATES: 500
NUMBER OF STEPS FOR CREATION RUN: 50000000


It is that OK when I set "begin 500 stop 500", but when I set "begin 50000500 stop 50000500" (I want to get the first frame of trajectory run2.dcd), I got the wrong pdb in which there is no coordinate information, just like:

ATOM 1 HT1 THR 1 ************************ 1.00 0.00 PROA
ATOM 2 HT2 THR 1 ************************ 1.00 0.00 PROA
ATOM 3 N THR 1 ************************ 1.00 1.00 PROA
ATOM 4 HT3 THR 1 ************************ 1.00 0.00 PROA
ATOM 5 CA THR 1 ************************ 1.00 1.00 PROA


Someone can help me to figure out this problem? Thanks


Last edited by zzlai; 05/09/11 10:09 PM.
Re: about reading frame from dcd one by one
zzlai #27445 05/09/11 10:48 PM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
I can think of at least 3 other ways to do this same task, all easier and more efficient (in terms of file I/O).

(1) process trajectory files one at a time; write out all coord sets as PDB, then process them, and then possibly delete the PDB files (I often use temporary storage, e.g. /tmp or /scratch).

(2) read all trajectory files in one pass; read coord set via TRAJ READ and write to PDB, followed by SYSTEM command to run the external program on each pass through the loop

(3) use "READ COOR FILE IFILE frame_number" to selectively read one coord set at a time; similar to original concept, but one can pass file and frame number as arguments, and use them to name the files; n == run number, j == frame number

charmm RUN:$n FRM:$j <getpdb.inp>getpdb.out
my codes n$n-$j.pdb


The CHARMM input would have

read coor file ifile @FRM name run@RUN.dcd
write coor pdb name n@RUN-@FRM.pdb


This is however, the least efficient, compared to (1) and (2) above.


Rick Venable
computational chemist

Re: about reading frame from dcd one by one
rmv #27446 05/09/11 11:40 PM
Joined: Mar 2010
Posts: 22
Z
zzlai Offline OP
Forum Member
OP Offline
Forum Member
Z
Joined: Mar 2010
Posts: 22
hi, thanks for the suggestions!
If I extract out all of the pdb files and analyse them, they will be extremely large. My computer doesn't have enough space to do that. My own codes are written by perl and fortran, I am not sure CHARMM can merge that. So I just simply use the third method. In the getpdb.inp, I just simply comment out

!TRAJECTORY FIRSTU 13 NUNIT 5 IREAD BEGIN 50000000
!TRAJ READ
!OPEN WRITE CARD UNIT 12 NAME temp123.pdb
!WRITE COOR PDB UNIT 12
!CLOSE UNIT 12

and add the scripts:

TRAJECTORY FIRSTU 13 NUNIT 5
read coor file ifile @FRM name run@RUN.dcd
write coor pdb name n@RUN_@FRM.pdb

and run like:

charmm RUN:2 FRM:2 <getpdb.inp>getpdb.out

Yet I got the following error:

CHARMM> TRAJECTORY FIRSTU 13 NUNIT 5
TRAJ: INITIATING READ OF A TRAJECTORY, OPTIONS;
FIRSTU = 13 NUNIT = 5 SKIP = 1

CHARMM> read coor file ifile @FRM name run@RUN.dcd
Parameter: FRM -> "2"
Parameter: RUN -> "2"
VOPEN> Attempting to open::run2.dcd::

***** LEVEL 0 WARNING FROM <MAINIO> *****
***** "READ...NAME...: OPEN" not possible.
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

Any suggestions are very welcome!

Re: about reading frame from dcd one by one
zzlai #27448 05/10/11 01:31 AM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
The TRAJ command isn't needed; either remove it or comment it out. It could be a problem if the file was already OPENed, i.e. the script also includes

open unit 13 file read name run@RUN.dcd

If so (I'm guessing), both that OPEN and the corresponding TRAJ command should be disabled.

Try the following changes--

open unit 11 read file name run@RUN.dcd
read coor file unit 11 ifile @FRM
close unit 11
write coor pdb name n@RUN-@FRM.pdb
* run @RUN fame @FRM
*


This uses the older 'READ COOR' syntax with an explicit UNIT, and adds a REMARK record to the PDB file.


Rick Venable
computational chemist

Re: about reading frame from dcd one by one
rmv #27453 05/10/11 05:18 AM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
Which CHARMM version is this?

Re: about reading frame from dcd one by one
Kenno #27454 05/10/11 07:24 AM
Joined: Sep 2003
Posts: 4,794
Likes: 2
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,794
Likes: 2
It is better to keep the discussion in one thread.
Are you sure that coor dmat cannot do what you want?


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: about reading frame from dcd one by one
rmv #27456 05/10/11 05:23 PM
Joined: Mar 2010
Posts: 22
Z
zzlai Offline OP
Forum Member
OP Offline
Forum Member
Z
Joined: Mar 2010
Posts: 22
hi, thanks for the suggestions!

At this time, I comment out

!open read unit 13 file name run1.dcd
!open read unit 14 file name run2.dcd
!open read unit 15 file name run3.dcd
!open read unit 16 file name run4.dcd
!open read unit 17 file name run5.dcd

!TRAJECTORY FIRSTU 13 NUNIT 5
!traj read
!read coor file ifile @FRM name run@RUN.dcd
!write coor pdb name n@RUN_@FRM.pdb

and added the scripts:

open unit 11 read file name run@RUN.dcd
read coor file unit 11 ifile @FRM
close unit 11
write coor pdb name n@RUN-@FRM.pdb
* run @RUN fame @FRM
*

and run like

charmm RUN:3 FRM:2 <getpdb.inp>getpdb.out

I got no error this time, but in the n3_2.pdb, I still cannot get the coordinate:

ATOM 1 HT1 THR 1 ************************ 1.00 0.00 PROA
ATOM 2 HT2 THR 1 ************************ 1.00 0.00 PROA
ATOM 3 N THR 1 ************************ 1.00 0.00 PROA
ATOM 4 HT3 THR 1 ************************ 1.00 0.00 PROA
ATOM 5 CA THR 1 ************************ 1.00 0.00 PROA
ATOM 6 CB THR 1 ************************ 1.00 0.00 PROA
ATOM 7 OG1 THR 1 ************************ 1.00 0.00 PROA
ATOM 8 HG1 THR 1 ************************ 1.00 0.00 PROA


My charmm version is c35b5

Any suggestions are very welcome!

Re: about reading frame from dcd one by one
zzlai #27458 05/10/11 07:11 PM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
The ***** usually implies numbers that are too large to be printed.

Do all the coord sets show this problem, or just some of them?

Try an alternate analysis, such as creating an ATOM XYZ time series for a single atom using all of the trajectory files.


Rick Venable
computational chemist

Re: about reading frame from dcd one by one
rmv #27459 05/10/11 08:12 PM
Joined: Mar 2010
Posts: 22
Z
zzlai Offline OP
Forum Member
OP Offline
Forum Member
Z
Joined: Mar 2010
Posts: 22
it seems like the coordinates are too large to be printed. When I check the pdb files which can be printed, it seems like the coordinates are too large, for example:

ATOM 1 HT1 THR 1 -882.680-611.079-891.411 1.00 0.00 PROA
ATOM 2 HT2 THR 1 -882.729-611.730-892.956 1.00 0.00 PROA
ATOM 3 N THR 1 -883.063-610.939-892.368 1.00 0.00 PROA
ATOM 4 HT3 THR 1 -884.094-610.884-892.244 1.00 0.00 PROA
ATOM 5 CA THR 1 -882.661-609.675-892.864 1.00 0.00 PROA
ATOM 6 CB THR 1 -883.484-608.541-892.379 1.00 0.00 PROA
ATOM 7 OG1 THR 1 -884.844-608.617-892.736 1.00 0.00 PROA
ATOM 8 HG1 THR 1 -885.202-607.891-892.220 1.00 0.00 PROA

it is very wierd, using this coordinate number to calculate RMSD, I still can get the reasonable vale, and when I put this pdb file into pymol, I find that the structure is still find. It doesn't crash. But my initial pdb is&#65306;

ATOM 1 HT1 THR 1 0.183 -8.780 3.376 1.00 0.00 PROA
ATOM 2 HT2 THR 1 -0.175 -7.127 3.221 1.00 0.00 PROA
ATOM 3 N THR 1 0.320 -7.933 2.788 1.00 1.00 PROA
ATOM 4 HT3 THR 1 1.336 -7.720 2.720 1.00 0.00 PROA
ATOM 5 CA THR 1 -0.215 -8.173 1.459 1.00 1.00 PROA
ATOM 6 CB THR 1 -1.679 -8.593 1.607 1.00 1.00 PROA
ATOM 7 OG1 THR 1 -1.668 -9.514 2.695 1.00 1.00 PROA
ATOM 8 HG1 THR 1 -2.560 -9.830 2.856 1.00 0.00 PROA

for each trajectory, at the end the simulation, I printed out a crd file. They are

final_1.crd

1 1 THR HT1 115.79449 -30.33708-111.85692 PROA 1 1.00000
2 1 THR HT2 116.90491 -29.08871-112.06959 PROA 1 1.00000
3 1 THR N 116.57857 -29.99956-112.45098 PROA 1 1.60000
4 1 THR HT3 117.29823 -30.71991-112.66258 PROA 1 1.00000
5 1 THR CA 116.05139 -29.98228-113.79038 PROA 1 2.36500
6 1 THR CB 117.01321 -29.04494-114.51789 PROA 1 2.36500
7 1 THR OG1 116.73869 -27.79922-113.84823 PROA 1 1.60000
8 1 THR HG1 117.37257 -27.17612-114.21090 PROA 1 1.00000


final_2.crd

1 1 THR HT1 -120.21912 24.98689 80.54568 PROA 1 1.00000
2 1 THR HT2 -120.63651 24.20858 82.04074 PROA 1 1.00000
3 1 THR N -120.97196 24.55674 81.11995 PROA 1 1.60000
4 1 THR HT3 -121.62629 25.32778 81.36275 PROA 1 1.00000
5 1 THR CA -121.51114 23.43157 80.43946 PROA 1 2.36500
6 1 THR CB -120.29825 22.62522 79.85330 PROA 1 2.36500
7 1 THR OG1 -119.56848 23.67411 79.19797 PROA 1 1.60000
8 1 THR HG1 -120.13022 23.95827 78.47319 PROA 1 1.00000

The first input file for the simulation is:

! Read topology and parameter files
open read card unit 10 name toph19.inp
read rtf card unit 10

open read card unit 20 name param19.inp
read para card unit 20

! Read PSF
open read unit 10 card name step1_pdbreader.psf
read psf unit 10 card

! Read COOR
open read unit 10 card name step1_pdbreader.pdb
read coor unit 10 pdb resid

!
! FACTS setup (param19)
!
! Protein dielectric constant
set diele 2.0
! Non-polar surface tension [cal/A^2]
set gamma 0.015

nbond nbxmod 5 atom cdiel eps @diele shift vatom vdistance vshift -
cutnb 9.0 ctofnb 7.5 ctonnb 6.5 e14fac 0.4 wmin 1.5

scalar wmain = radius
scalar wmain set 1.0 selection (type h*) end

facts tcps 19 teps @diele gamm @gamma tavw


ENERGY

!
! short minimization and Langevin dynamics
! NOTE: CHANGE ISEED EACH TIME FOR RESTART IN LANGEVIN DYNAMICS
!

ENERGY
! short minimization and Langevin dynamics
! NOTE: CHANGE ISEED EACH TIME FOR RESTART IN LANGEVIN DYNAMICS
!

mini sd nstep 50 nprint 10 step 0.005 inbfrq -1
mini abnr nstep 50 nprint 10 step 0.005 inbfrq -1

set fbeta = 0.15
set temp = 330.0

scalar fbeta set @fbeta
scalar fbeta set 0.0 select type H* end

SHAKE BONH PARAm TOL 1.0e-6

open write unit 12 card name run1.rst
open write unit 13 file name run1.dcd

DYNAMICS LANGEVIN start nstep 12500000 timestp 0.002 -
nprint 1000 iprfrq 1000 isvfrq 1000 ntrfrq 500 -
inbfrq -1 imgfrq -1 ihbfrq 0 ilbfrq 0 -
firstt @temp finalt @temp iseed 8304 -
tbath @temp rbuf 0.0 -
iunread -1 iunwrite 12 iuncrd 13 iunvelo -1 -
nsavcrd 500 nsavvelo 0 -
iasvel 1



open write card unit 10 name final_1.crd
write coor card unit 10

stop


The second (restart)input file is:


! Read topology and parameter files
open read card unit 10 name toph19.inp
read rtf card unit 10

open read card unit 20 name param19.inp
read para card unit 20



! Read PSF
open read unit 10 card name step1_pdbreader.psf
read psf unit 10 card

! Read COOR
open read unit 10 card name final_1.crd
read coor unit 10 card

!
! FACTS setup (param22)
!
! Protein dielectric constant
set diele 2.0
! Non-polar surface tension [cal/A^2]
set gamma 0.015

nbond nbxmod 5 atom cdiel eps @diele shift vatom vdistance vshift -
cutnb 9.0 ctofnb 7.5 ctonnb 6.5 e14fac 0.4 wmin 1.5

scalar wmain = radius
scalar wmain set 1.0 selection (type h*) end

facts tcps 19 teps @diele gamm @gamma tavw


ENERGY

!
! short minimization and Langevin dynamics
! NOTE: CHANGE ISEED EACH TIME FOR RESTART IN LANGEVIN DYNAMICS
!

!mini sd nstep 50 nprint 10 step 0.005 inbfrq -1
!mini abnr nstep 50 nprint 10 step 0.005 inbfrq -1

set fbeta = 0.15
set temp = 330.0

scalar fbeta set @fbeta
scalar fbeta set 0.0 select type H* end

SHAKE BONH PARAm TOL 1.0e-6

open read unit 11 card name run1.rst
open write unit 12 card name run2.rst
open write unit 13 file name run2.dcd

DYNAMICS LANGEVIN restart nstep 12500000 timestp 0.002 -
nprint 1000 iprfrq 1000 isvfrq 0 ntrfrq 500 -
inbfrq -1 imgfrq -1 ihbfrq 0 ilbfrq 0 -
firstt @temp finalt @temp -
tbath @temp rbuf 0.0 -
iunread 11 iunwrite 12 iuncrd 13 iunvelo -1 -
nsavcrd 500 nsavvelo 0 -
iasvel 1


ENERGY

open write card unit 10 name final_2.crd
write coor card unit 10

stop


Something wrong?

Last edited by zzlai; 05/10/11 08:17 PM.
Re: about reading frame from dcd one by one
zzlai #27461 05/10/11 09:26 PM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
The protein is simply translating in space. Periodic boundary conditions could be used to keep the coordinates within a specific volume of space, avoiding the problem of coordinates which are too large for the PDB format.


Rick Venable
computational chemist

Page 1 of 2 1 2

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.4
(Release build 20200307)
Responsive Width:

PHP: 5.6.33-0+deb8u1 Page Time: 0.012s Queries: 35 (0.005s) Memory: 0.9968 MB (Peak: 1.1298 MB) Data Comp: Off Server Time: 2020-09-30 22:33:07 UTC
Valid HTML 5 and Valid CSS