CHARMM Development Project Forums CHARMM Community Script Archive How to obtain the atomicforce on each protein atom

Hi all,

I am a beginner on charmm-script and I am recently trying to obtain the atomic force acting on each protein atoms in my simulation. I read from CHARMM docs that I can calculate the force wity 'ENERgy FORCE' command but I am not sure how I can get the forces, like, output it to a single file.

Here is my current charmm script for reading in dcd file and loop over trajectories:




! Open and read DCD file
open read unit 90 file name dyn.dcd

open write unit 14 form name potential_energy.dat

! Open file to write forces
!open write unit 15 form name atomic_forces.dat

set i 1
TRAJECTORY IREAD 90 BEGIN 40000 SKIP 40000
label loop
Traj read
ENERgy FORCE
write title unit 14
* Frame: @i Energy: ?ener Box_Size: @A @B @C FORCE: ?F
*
! ! Print forces on all atoms
! set n 1
! label forces_loop
! write title unit 14
! * Atom: @n Force: ?FARRAY(@i,1) ?FARRAY(@i,2) ?FARRAY(@i,3)
! *
! incr n by 1
! if n le ?NATOM goto forces_loop


incr i by 1
!Check if there are more frames in trajectory
if i lt 5000 goto loop



STOP
!Close file unit for writing energies
close unit 14
close unit 90

You do not state what goes wrong with your attempt. Do you get any error messages?

Start with just one coordinate set (ie, just a coordinate file, not a trajectory) and see to it that you can get the information you need for that. Only then add the loop over the trajectory.
I am not sure where you got all your commands from, some seem to be wild guesses (does not work).


Basically
READ RTF ...
READ PARAM...
READ PSF ....
any special setup, e.g. PBC, restraints, ....

READ COOR CARD NAME MY_SYSTEM.CRD
ENERGY ...
COOR FORCE COMP ! copy forces to the comparison coordinates
COOR DIST WEIGH COMP ! put magnitudes in the weighting array
WRITE COOR CARD COMP NAME FORCES_MY_SYSTEM. DAT
* force on each atom. magnitude is in weighting array. CHARMM variables can be added here
*

This produces quite a bit of output; most of it is useful since it retains the residue and atom names. You can trim it down in a couple of ways.
1/ write out the forces for only a (small) subset of atoms in the system - add an atom selection to the write command. Note that you probably want the whole system to be present for the energy/force evaluation.
2/ extract just the magnitude of the force on each atom (requires some manipulations using SCALAr commands). I am not sure how easy it would be to get the id of each atom here.

Hi Lennart,

Sorry I did not make myself clear, so in my current script, at the line '* Frame: @i Energy: ?ener Box_Size: @A @B @C FORCE: ?F', I can use ?ener to visit the potential energy calculated from my system and save them in the a file and in similar way I want to save the atomic force on each protein atoms, so there is a 'FORCE: ?F' term. I am not sure what I should but instead of ?F cause ?F will not give me the atomic force

I have given you the commands you need. Try them.
There is no ?F or ?FARRAY internal variable, and no ENERgy FORCe command (this is what I mean when I say that you cannot invent command and options yourself). Take the error message seriously!

Thank you so much! I tried it and it works well. But I think this command copies the forces on all atoms in the system, including water, ions and protein atoms. Would you mind give me some more instructions on how to print only protein atomic forces? I read your pre-post that I can modify the write command to select the atom groups, so do I modify it in the line of WRITE command? Or do I have to define a group of atoms before the force calculation part?

Thanks!

Well, depends on what you want - I cannot do the science for you. Check with your supervisor.
The force on an atom during the simulation comes from interactions with all the rest of the system (including PBC images), so if this is the force you are interested in you have to have the same setup of your system and the same options for the energy evaluation as you had during the simulation.
Another issue is for which atoms you want to get the forces.
The WRITE COOR command takes an atom selection (io.info, WRITe section), e.g. WRITE COOR CARD NAME MYFILE.CRD SELE SEGI PROT END
subst.info lists the internal variables that you can access using ?name. The CRYSTAL section has information about lattice.

Hi, sorry to bother again.

I tried the command you gave me before and I got the force for one frame printed out.

I also tried another way in which I set the parameter NSAVX in the DYNAMIC part of CHARMM input file, which allows me to save the force as the simulation going. However, the force I saved in this way is different from the force I got using your command, do you know why this happens?

Following is the DYNA part of charmm input and the command I use to print force from reading traj:

CHARMM INP:
dynamics cpt leap start timestep 0.001 nstep @ndyn -
firstt @temperature finalt @temperature twindl -10.0 twindh 10.0 -
iasvel 5 iasors 5 iunr -8 iunwri 12 -
nsavc 4000 nsavv 4000 nsavx 4000 MXYZ 3 nprint 1000 ntrfrq 1000 -
pconstant pmass 500.0 pref 1.0 pgamma 20.0 -
hoover reft @temperature tmass 1000.0 tbath @temperature -
iuncrd 13 iunvel 13 iunxyz 20 isvfrq @resfrq -
echeck -1.0 ! iseed 363520106 363520106 363520106 363520106

Command you gave me:

TRAJECTORY IREAD 90 BEGIN 0 SKIP 1
label loop
Traj read
ENERgy
COOR STAT
COOR FORCE COMP
COOR DIST WEIGH COMP ! put magnitudes in the weighting array
WRITE COOR CARD COMP NAME FORCES_MY_SYSTEM_@i.DAT

Check that ALL options relevant for the energy evaluation are the same in both cases. Cutoffs, PBC, restraints, ...
The information you have posted is not sufficient for us to tell if this is the case.
You use the same I/O-unit for velocities and coordinates, is this intentional? You should not normally have to set ISVFRQ. It is possible, but messy, to restart a simulation from a restart file that was saved during the run.
Why do you turn off the energy check (echeck -1)?

I think they are the same, I mean for ALL options relevant for the energy evaluation are the same in both cases. Because the potential energy and box size for both cases are the same for each frame.

If the energy is the same the forces should also be the same. How different are they?