Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
Calculate the correlation between the original B-f
#24000 04/14/10 02:19 AM
Joined: Oct 2009
Posts: 22
J
Jun Offline OP
Forum Member
OP Offline
Forum Member
J
Joined: Oct 2009
Posts: 22
Hi, all.

I am JunKoo and am quite new here.
I want to calculate the correlation between B-factor of original protein data(PDB file) and that of calculated data using normal modes.
So, I want to calculate rmsf for each atom using normal modes.

I used the vibran command and diag worked successfully.
and I am having hard time using rmsf-residue-str.

My questions is
1) Do I need to use rmsf-residue.str ? Is there any other way?

2) If I use rmsf-residue-str, how to choose and make trajectory in order to get traj.cor ?
Does rmsf depend on the trajectory that I choose ?

3) When I use dynamics.inp, I get the traj.cor
When I run charmm, I got the following error.

[jun@protein sample]$ charmm <rmsf1not.inp> rmsf1not.out
sue : unformatted io not allowed
apparent state: unit 90 named 1not.psf
lately reading sequential unformatted external IO
Aborted





Thank you.



Last edited by Jun; 04/14/10 02:27 AM.
Re: Calculate the correlation between the original B-f
Jun #24002 04/14/10 03:20 AM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
The rmsf-residue.str script is designed for analyzing RMSF from a trajectory via COOR DYNA, but may also use fluctuations from normal modes. However, the FLUCT ATOM command in VIBRAN provides the fluctuation vector in the x,y,z coord positions, while COOR DYNA puts the isotropic fluctuations into the weighting array, so a few extra steps may be needed. A couple options are: [1] use SCALAR commands to compute isotropic fluctuations, or [2] write a normal mode trajectory and process it with COOR DYNA.

The rmsf-residue.str script also seems to have an error for the PSF reading example; the CARD keyword is omitted.


Rick Venable
computational chemist

Re: Calculate the correlation between the original B-f
rmv #24008 04/14/10 03:51 PM
Joined: Oct 2009
Posts: 22
J
Jun Offline OP
Forum Member
OP Offline
Forum Member
J
Joined: Oct 2009
Posts: 22
Thank you for the quick reply.


When I use the FLUCT ATOM command, I got the following result.

VIBRAN> fluctuations atom

NORMAL MODE FLUCTUATIONS

ATOM FLUCTUATIONS WILL BE COMPUTED
NO MASS WEIGHTING WILL BE USED
MODE 1 FREQ= -15.15 FACTOR= 1.000000
RMS= 0.000565 TOTAL= 0.000565 %INCR= 100.000000
MODE 2 FREQ= 0.00 FACTOR= 1.000000
RMS= 0.000348 TOTAL= 0.000913 %INCR= 38.079117
MODE 3 FREQ= 0.00 FACTOR= 1.000000
RMS= 0.000348 TOTAL= 0.001260 %INCR= 27.577752
MODE 4 FREQ= 0.00 FACTOR= 1.000000
RMS= 0.000348 TOTAL= 0.001608 %INCR= 21.616427
MODE 5 FREQ= 4.85 FACTOR= 1.000000
RMS= 0.000384 TOTAL= 0.001992 %INCR= 19.271940
MODE 6 FREQ= 5.55 FACTOR= 1.000000
RMS= 0.000378 TOTAL= 0.002370 %INCR= 15.967969
MODE 7 FREQ= 7.03 FACTOR= 1.000000
RMS= 0.000417 TOTAL= 0.002787 %INCR= 14.964332
MODE 8 FREQ= 11.52 FACTOR= 1.000000
RMS= 0.000400 TOTAL= 0.003187 %INCR= 12.545851
MODE 9 FREQ= 12.87 FACTOR= 1.000000
RMS= 0.000448 TOTAL= 0.003635 %INCR= 12.332399
MODE 10 FREQ= 16.37 FACTOR= 1.000000
RMS= 0.000425 TOTAL= 0.004060 %INCR= 10.464939


.............
.............
.............


MODE 525 FREQ= 3514.85 FACTOR= 1.000000
RMS= 0.002632 TOTAL= 0.721874 %INCR= 0.364572
MODE 526 FREQ= 3536.48 FACTOR= 1.000000
RMS= 0.002542 TOTAL= 0.724416 %INCR= 0.350864
MODE 527 FREQ= 3647.68 FACTOR= 1.000000
RMS= 0.002651 TOTAL= 0.727067 %INCR= 0.364588
MODE 528 FREQ= 3679.18 FACTOR= 1.000000
RMS= 0.002657 TOTAL= 0.729724 %INCR= 0.364158
AVERAGE FLUCTUATION OVER SELECTED ATOMS IS 0.8542

VIBRAN> end


My question is
1) To calculate the correlation, I think I need to have the same number of atoms between pdb files and my results.
In PDB file, the number of atom = 96
In PSF file, the number of atom = 176
In the result, I only have the number of modes = 528

So, Do the modes 1,2,3 correspond to atom 1 ?
then I need to average the RMS value of mode 1,2,3 for atom1.
Is it correct ?

2) Because the number of atoms in PDB file and the number of atoms in PSF file differs, how to compare the B-factor of atoms in PDB file and the B-factor of atoms in PSF file?


Thank you in advance.

Re: Calculate the correlation between the original B-f
Jun #24009 04/14/10 06:48 PM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
Modes 1-6 are translation and rotation, and should all be zero for a structure which has been adequately minimized in preparation for a normal mode calculation. You have some more work to do.

I'm guessing that the PDB file may lack H atoms (you should verify this), so you should probably limit the comparison to the heavy atoms.


Rick Venable
computational chemist

Re: Calculate the correlation between the original B-f
rmv #24015 04/15/10 02:49 AM
Joined: Oct 2009
Posts: 22
J
Jun Offline OP
Forum Member
OP Offline
Forum Member
J
Joined: Oct 2009
Posts: 22
I was trying to minimize it, but at some point the frequency for rotation does not go to zero.(frequency for translation is zero)

my input for minimization

---------------------------------------------
read rtf card name ~/charmm/toppar/top_all22_prot.inp
read param card name ~/charmm/toppar/par_all22_prot.inp

! Get psf and coordinates
read psf card name 1not.psf
read coor card name 1not.crd

coor copy comp

! Harmonic constraints on all protein atoms
cons harm force 20.0 sele segid 1not end

! First used Steepest Descent (a gentle minimizer)
minimize sd nstep 1000 cdie eps 1.0 fshift vshift cutnb 13.0 ctofnb 12.0

! How different is the structure now?
coor orie rms sele segid 1not end
! And if we look at the waters instead?
coor orie rms sele segid xwat end

! Reduce harmonic constraints and minimize some more
cons harm force 10.0 sele segid 1not end
minimize abnr nstep 1000
coor orie rms sele segid 1not end
! turn off harmonic constraints
cons harm force 0.0 sele segid 1not end
! then small constraints on backbone only
cons harm force 5.0 sele segid 1not .and. -
(type C .or. type N .or. type CA .or. type O) end
minimize abnr nstep 2000
coor orie rms sele segid 1not end
cons harm force 0.0 sele segid 1not end
minimize abnr nstep 3000
coor orie rms sele segid 1not end

write coor card name 1not_min.crd

----------------------------------------------------
output


SELRPN> 176 atoms have been selected out of 176
CENTER OF ATOMS BEFORE TRANSLATION -1.80604 -0.71568 -7.27346
CENTER OF REFERENCE COORDINATE SET -1.80604 -0.71568 -7.27346
NET TRANSLATION OF ROTATED ATOMS 0.00000 0.00000 0.00000
ROTATION MATRIX
0.999998 0.000284 0.002027
-0.000282 1.000000 -0.000741
-0.002027 0.000741 0.999998
AXIS OF ROTATION IS -0.340445 -0.931242 0.129944 ANGLE IS 0.12

TOTAL SQUARE DIFF IS 281.1198 DENOMINATOR IS 176.0000
THUS RMS DIFF IS 1.263832
ALL COORDINATES ORIENTED IN THE MAIN SET BASED ON SELECTED ATOMS.


------------------------------------------------
my output for vibran is

VIBRAN> fluctuations atom

NORMAL MODE FLUCTUATIONS

ATOM FLUCTUATIONS WILL BE COMPUTED
NO MASS WEIGHTING WILL BE USED
MODE 1 FREQ= 0.00 FACTOR= 1.000000
RMS= 0.000348 TOTAL= 0.000348 %INCR= 100.000000
MODE 2 FREQ= 0.00 FACTOR= 1.000000
RMS= 0.000348 TOTAL= 0.000695 %INCR= 50.000000
MODE 3 FREQ= 0.00 FACTOR= 1.000000
RMS= 0.000348 TOTAL= 0.001043 %INCR= 33.333333
MODE 4 FREQ= 4.00 FACTOR= 1.000000
RMS= 0.000383 TOTAL= 0.001426 %INCR= 26.872102
MODE 5 FREQ= 5.04 FACTOR= 1.000000
RMS= 0.000368 TOTAL= 0.001794 %INCR= 20.506263
MODE 6 FREQ= 6.12 FACTOR= 1.000000
RMS= 0.000398 TOTAL= 0.002191 %INCR= 18.147414



1) Would you give me some suggestions ?

2) When I choose the heavy atoms, I used the following command.

cons fix sele prop mass .ge. 12.012 end

When I use the cutoff 12.012, I have 42 atoms.
When I use the cutoff 12.011, I have 97 atoms.

How can I choose 96 atoms? Because I have 96atoms in the PDB file.


Thank you.

Re: Calculate the correlation between the original B-f
Jun #24018 04/15/10 04:24 AM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
I don't understand what the CONS FIX command is for.

Limit the SD minimization to 20-50 steps, just enough to handle any bad atom contacts; monitor the VDW energy term for this.

Heavy atoms are non-hydrogen atoms (sele .not. hydrogen end)

CHARMM may have added a missing heavy atom, such as one of the C-term carboxylate O atoms.


Rick Venable
computational chemist

Re: Calculate the correlation between the original B-f
rmv #24022 04/15/10 04:56 PM
Joined: Oct 2009
Posts: 22
J
Jun Offline OP
Forum Member
OP Offline
Forum Member
J
Joined: Oct 2009
Posts: 22
I used the 20 steps for SD. I have the following output.


NSTEP = 20 NPRINT = 10
STEP = 0.0200000 TOLFUN = 0.0000000
TOLGRD = 0.0000000 TOLSTP = 0.0000000

MINI MIN: Cycle ENERgy Delta-E GRMS Step-size
MINI INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
MINI EXTERN: VDWaals ELEC HBONds ASP USER
---------- --------- --------- --------- --------- ---------
MINI> 0 -14.35464 0.00000 9.05141 0.02000
MINI INTERN> 12.48736 30.89169 2.77434 64.60245 0.66050
MINI EXTERN> -1.52277 -124.24820 0.00000 0.00000 0.00000
---------- --------- --------- --------- --------- ---------
MINI> 10 -44.45200 30.09735 0.91057 0.00156
MINI INTERN> 6.84458 23.99652 2.14278 63.64631 1.28447
MINI EXTERN> -15.33326 -130.60277 0.00000 0.00000 0.00000
MINI CONSTR> 3.56937 0.00000 0.00000 0.00000 0.00000
---------- --------- --------- --------- --------- ---------
MINI> 20 -46.54647 2.09447 0.89600 0.00070
MINI INTERN> 6.88212 24.12926 2.17575 63.50708 1.21484
MINI EXTERN> -17.84230 -132.67036 0.00000 0.00000 0.00000
MINI CONSTR> 6.05715 0.00000 0.00000 0.00000 0.00000
---------- --------- --------- --------- --------- ---------

STEEPD> Minimization exiting with number of steps limit ( 20) exceeded.

STPD MIN: Cycle ENERgy Delta-E GRMS Step-size
STPD INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
STPD EXTERN: VDWaals ELEC HBONds ASP USER
STPD CONSTR: HARMonic CDIHedral CIC RESDistance NOE
---------- --------- --------- --------- --------- ---------
STPD> 20 -46.54647 2.09447 0.89600 0.00084
STPD INTERN> 6.88212 24.12926 2.17575 63.50708 1.21484
STPD EXTERN> -17.84230 -132.67036 0.00000 0.00000 0.00000
STPD CONSTR> 6.05715 0.00000 0.00000 0.00000 0.00000

1) It seems like that the VDW is -17 at the end. Do I have any bad atoms?


2) I like to write non-hydrogen atoms and their RMSF values in external file. What kind of command do I use and where do I put the command? (Inside of vibran command or outside of vibran command ) I read IO.doc, however, I couldn't understand.

This is my input for selecting non-hydrogen atoms.

vibran nmod 528
diag
fluctuations atom
!sele .not. hydrogen end
!coor orie rms sele segid 1not end
!write coor card name 1not.dat
end

!open unit 12 write coor name 1not.dat sele .not. hydrogen end

!coor copy sele .not. hydrogen end



3) It seems like CHARMM may add a missing heavy atom because when I choose non-hydrogen atoms, the number of atoms(=97) is bigger than the number of atoms(=96) in PDB file. How can I exclude the atom like (C-term carboxylate O atoms) ? I was seeing the topology file to figure out.



Thank you for your help.

Re: Calculate the correlation between the original B-f
Jun #24024 04/15/10 06:37 PM
Joined: Sep 2003
Posts: 8,499
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,499
1) Atoms which are initially too close to each other can lead to very high VDW energy, but that does seem to be the case here.

2) As noted above, to get isotropic RMSF you must either: [1] use SCALAR commands (scalar.doc) to compute isotropic fluctuations, or [2] write a normal mode trajectory (WRITE command in VIBRAN) and process it with COOR DYNA.

3) Try comparing your CHARMM coordinate file to the original PDB file, or (after import) maybe write out a heavy atom PDB file via

write coor pdb name hvy1not.pdb sele .not. hydrogen end


Rick Venable
computational chemist

Re: Calculate the correlation between the original B-f
rmv #24027 04/15/10 07:56 PM
Joined: Oct 2009
Posts: 22
J
Jun Offline OP
Forum Member
OP Offline
Forum Member
J
Joined: Oct 2009
Posts: 22
Thank you very much.

For 2), which method do you recommend?
To use rmsf.residue.str, I need to create 2 trajectories.
I made trajectory using dynamic.inp Is it OK? or
I need to use write command in VIBRAN?

When I made trajectory, I have the following output.

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

TITLE> *FILENAME: DYNAMICS.INP
TITLE> *PURPOSE: VERY BASIC MD EXAMPLE OF PROTEIN WITH XTAL WATERS IN VACUUM
TITLE> *AUTHOR: LENNART NILSSON, KAROLINSKA INSTITUTET, OCTOBER 2003
TITLE> * DATE: 4/14/10 1:26:56 CREATED BY USER: jun
TITLE> *

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

***** ERROR ***** FILES DO NOT MERGE
UNIT 52 BEGINS AT 100 INSTEAD OF 10100
*** LEVEL -1 WARNING *** BOMLEV IS 0
BOMLEV HAS BEEN SATISFIED. TERMINATING.


Q. If I use write command in VIBRAN, would you give me some examples? After creating the trajectory, I need to use rmsf-residue.str, correct?


Thank you always.

Re: Calculate the correlation between the original B-f
Jun #24028 04/15/10 08:14 PM
Joined: Sep 2003
Posts: 4,794
Likes: 2
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,794
Likes: 2
Your post is fairly incoherent; this is not good, neither for your science, nor for our ability to help you

I think you need to clarify to yourself what it is you want to achieve, and then consider carefully how to do this using CHARMM commands. For example, you seem to be a bit confused about different numbers of atoms in the PDB file and in the corresponding CHARMM structure; this is not really a CHARMM issue, but rather an issue of knowing what you are doing. You need to be able to resolve such problems before asking specific questions as to how to implement your intentions in a CHARMM procedure.

Note that there is no "rmsf.residue.str" example in the Script Archive (attention to detail is extremely important when dealing with computers!).

It is OK to use a file called dynamic.inp to create your trajectory - the name of the file does not really matter. Unfortunately the contents of the file are very important.

rmsf-residue.str does not use two trajectories. Why are you so focused on rmsf-residue.str? It has very little (nothing!) to do with correlating B-factors from normal mode analysis and Xr-ray protein crystallography.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
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: 34 (0.007s) Memory: 0.9971 MB (Peak: 1.1269 MB) Data Comp: Off Server Time: 2020-09-30 20:52:31 UTC
Valid HTML 5 and Valid CSS