CHARMM Development Project Forums User Discussion & Questions Setup, I/O, and Basic questions error Computi RMSD and Rg using CORREL

Hi CHARMM users,

I while runing an analisys script for rmsd-rgyr-trajectory correl. I had the following problem:

Please help me to find below the problematic part of the input file and below this the CHARMM output.

Thanks for your help


input
* Compute RMSD and Rg for a MD trajectory: using CORREL tools
*

!Read topology and parameter (aka force field)
read rtf card name /home/lherrera/proyectos/namd/top_all22_prot_cmap_cu.rtf
read para card name /home/lherrera/proyectos/namd/par_all22_prot_cmap_cu.prm

!Read PSF and initial/reference coordinates
read psf card name /home/lherrera/proyectos/charmm/charmm-analysis/files/3gb1_solv.psf
read coor pdb resid name /home/lherrera/proyectos/charmm/charmm-analysis/files/3gb1_solv.pdb
coor copy comp

!Open file unit for trajectory
open read unit 13 file name /home/lherrera/proyectos/charmm/charmm-analysis/files/nptprod.pdb
traj query unit 13

!Open the output file and write header
open write unit 11 card name /home/lherrera/proyectos/charmm/charmm-analysis/files/outleo/rmsd-rgyr-correl.dat
write title unit 11
*#time rmsd rgyr (processed by CORREL)
*

!Invoke CORREL mode
correl maxtime ?nfile

!Now we are in CORREL mode and have access to CORREL subcommands

!request RMS with orient and radius of gyration
enter v1 rms orient
enter v2 gyra

!specify the trajectory and atom selection to process
traj firstu 13 nunit 1 select .not. resn tip3 end

!combine v1 and v2 to allow both time series to be output to one file
edit v1 veccod 2

!write the time series to a file
write v1 dumb time unit 11

! exit CORREL mode
end

stop



output
....
....
....
CHARMM> write title unit 11
RDTITL> *#TIME RMSD RGYR (PROCESSED BY CORREL)
RDTITL> *

CHARMM>

CHARMM> !Invoke CORREL mode
CHARMM> correl maxtime ?nfile
RDCMND substituted energy or value "?NFILE" to "1095649606"

NONBOND OPTION FLAGS:
ELEC VDW ATOMs CDIElec SHIFt VATOm VSWItch
BYGRoup NOEXtnd NOEWald
CUTNB = 14.000 CTEXNB =999.000 CTONNB = 10.000 CTOFNB = 12.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.
<MAKINB> with mode 5 found 18645 exclusions and 2239 interactions(1-4)
<MAKGRP> found 762 group exclusions.
Generating nonbond list with Exclusion mode = 5
== PRIMARY == SPACE FOR 4910265 ATOM PAIRS AND 0 GROUP PAIRS
== PRIMARY == SPACE FOR 7365417 ATOM PAIRS AND 0 GROUP PAIRS

General atom nonbond list generation found:
7313984 ATOM PAIRS WERE FOUND FOR ATOM LIST
404382 GROUP PAIRS REQUIRED ATOM SEARCHES


***** LEVEL -4 WARNING FROM <ALLHP> *****
***** Bad length was passed to FREEHP
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5



/---------\
/ \
/ \
/ \
! XXXX XXXX !
! XXXX XXXX !
! XXX XXX !
! X !
--\ XXX /--
! ! XXX ! !
! ! ! !
! I I I I I !
! I I I I !
\ /
-- --
\---/
XXX XXX
XXXX XXXX
XXXXX XXXXX
XXX XXX
XXX XXX
XXXXX
XXX XXX
XXX XXX
XXX XXX
XXXXX XXXXX
XXXX XXXX
XXX XXX


Execution terminated due to the detection of a fatal error.

ABNORMAL TERMINATION
MAXIMUM STACK SPACE USED IS 307018
STACK CURRENTLY IN USE IS 62
MOST SEVERE WARNING WAS AT LEVEL -4
HEAP PRINTOUT- HEAP SIZE 10240000
SPACE CURRENTLY IN USE IS 7733386
MAXIMUM SPACE USED IS 7733386
FREE LIST
PRINHP> ADDRESS: 1 LENGTH: 2506082 NEXT: 10239469
PRINHP> ADDRESS: 10239469 LENGTH: 532 NEXT: 0

$$$$$ JOB ACCOUNTING INFORMATION $$$$$
ELAPSED TIME: 2.91 SECONDS
CPU TIME: 2.91 SECONDS

Unless the binary trajectory file has 1095649606 frames in it, there may be a problem with that file; check the TRAJ QUERY output more carefully. Was the trajectory file produced with CHARMM?

my dcd file has 101 frames, how can I solve them?
thanks.

A few notes and suggestions--

• the file "nptprod.pdb" has an odd extension for a binary trajectory file; usually .trj or .dcd are used
• when not computing energy, use of the NOUPDATE keyword is highly recommended as part of the initial CORREL command
• specify the MAXTIME value explicitly, to e.g. 200

CORREL MAXT 200 NOUPDATE

You may still have trouble, esp. if NAMD is involved; in general, unless the initial PSF was created with CHARMM and exported as an XPLOR PSF for use with NAMD, you can expect problems.

There can also be binary file incompatibilities, if the trajectory was produced on a machine with a different architecture (big- vs. little-endian). CHARMM itself does not handle that very well at all.

hello to seem it again...
I solved problem changing 1095649606 by 101 NOUPDATE (I number of frames) is correct this? , but now I obtain another error, rotten to help me to solve this new misfortune.
thanks in advance


CHARMM> !Invoke CORREL mode
CHARMM> correl maxtime 101 NOUPDATE

CORREL>

CORREL> !Now we are in CORREL mode and have access to CORREL subcommands

CORREL>

CORREL> !request RMS with orient and radius of gyration

CORREL> enter v1 rms orient

CORREL> enter v2 gyra

CORREL>

CORREL> !specify the trajectory and atom selection to process

CORREL> traj firstu 13 nunit 1 select .not. resn tip3 end
SELRPN> 855 atoms have been selected out of 17088
The following time series will be filled:
V1
V2


***** ERROR ***** HEADERS DO NOT MATCH RK COOR CORD
*** LEVEL -1 WARNING *** BOMLEV IS 0
BOMLEV HAS BEEN SATISFIED. TERMINATING.



/---------\
/ \
/ \
/ \
! XXXX XXXX !
! XXXX XXXX !
! XXX XXX !
! X !
--\ XXX /--
! ! XXX ! !
! ! ! !
! I I I I I !
! I I I I !
\ /
-- --
\---/
XXX XXX
XXXX XXXX
XXXXX XXXXX
XXX XXX
XXX XXX
XXXXX
XXX XXX
XXX XXX
XXX XXX
XXXXX XXXXX
XXXX XXXX
XXX XXX


Execution terminated due to the detection of a fatal error.

ABNORMAL TERMINATION
MAXIMUM STACK SPACE USED IS 119624
STACK CURRENTLY IN USE IS 68420
MOST SEVERE WARNING WAS AT LEVEL -1
HEAP PRINTOUT- HEAP SIZE 10240000
SPACE CURRENTLY IN USE IS 97182
MAXIMUM SPACE USED IS 97182
FREE LIST
PRINHP> ADDRESS: 1 LENGTH: 10142818 NEXT: 0

$$$$$ JOB ACCOUNTING INFORMATION $$$$$
ELAPSED TIME: 1.97 SECONDS
CPU TIME: 1.82 SECONDS

CHARMM does not recognize your file as a valid binary trajectory; see my above response.

excuse me, I was confused my file dcd.bin with dcd normal ….
how can I make the same calculation for two residues (bridge salts)?
what is the sentence that I would have to include in my script? thanks in advance.

You want to do the same calculation (RMSD and Rg vs time) for two residues forming a salt bridge? You can select the two residues of interest in the atom selection (select.doc) on the TRAJ command in the correl module (correl.doc).

how ?
I had the other error:
CHARMM> !Open file unit for trajectory
CHARMM> open read unit 13 file name /home/lherrera/proyectos/charmm/charmm-analysis/files/nptprod.dcd
VOPEN> Attempting to open::/home/lherrera/proyectos/charmm/charmm-analysis/files/nptprod.dcd::
OPNLGU> Unit 13 opened for READONLY access to /home/lherrera/proyectos/charmm/charmm-analysis/files/nptprod.dcd

CHARMM> traj query unit 13

READING TRAJECTORY FROM UNIT 13
NUMBER OF COORDINATE SETS IN FILE: 100
NUMBER OF PREVIOUS DYNAMICS STEPS: 50
FREQUENCY FOR SAVING COORDINATES: 50
NUMBER OF STEPS FOR CREATION RUN: 5000

NUMBER OF DEGREES OF FREEDOM: 35028
NUMBER OF FIXED ATOMS: 0
THE INTEGRATION TIME STEP (PS): 0.0020
THE FILE CONTAINS CRYSTAL UNIT CELL DATA
THE FILE DOES NOT CONTAIN 4-D DATA

CHARMM>

CHARMM> !Open the output file and write header
CHARMM> open write unit 11 card name /home/lherrera/proyectos/charmm/charmm-analysis/files/outleo/rmsd2redrgcor.dat
VOPEN> Attempting to open::/home/lherrera/proyectos/charmm/charmm-analysis/files/outleo/rmsd2redrgcor.dat::
OPNLGU> Unit 11 opened for WRITE access to /home/lherrera/proyectos/charmm/charmm-analysis/files/outleo/rmsd2redrgcor.dat

CHARMM> write title unit 11
RDTITL> *#TIME RMSD RGYR (PROCESSED BY CORREL)
RDTITL> *

CHARMM>

CHARMM> sele resid 1 .and resid 10

***** LEVEL 0 WARNING FROM <CHARMM> *****
***** Unrecognized command: SELE
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

Thanks

I don't think that Rg for just a couple residues is meaningful.

Major syntax errors here; you have some reading to do. Within the CORREL module (correl.doc), a few of the time series types use or allow an atom selection (select.doc).

traj firstu 13 nunit 1 select .not. resn tip3 end

The atom selection in the above command within CORREL from your script is meant to select the entire protein, assuming that only water (TIP3) is present. You need to change that and get the syntax and the logic correct, e.g.

traj firstu 13 nunit 1 select segid A .and. ( resid 1 .or. resid 10 ) end

I've assumed the protein is SEGID A; each segment can have a RESID 1, so the RESID alone is not unique for more than a single segment,