Page 1 of 2 1 2 >
Topic Options
#7796 - 08/10/05 12:12 PM Is this possible in CHARMM (adding a routine) ?
cola Offline
Forum Member

Registered: 07/28/04
Posts: 113
Dear All,

Let's assume that I have three atoms.
I'd like to get each position and
calculate the normal vector to a plane which is made by these three atoms.
(The normal vector will be changed during the simulation.)
Then, I'd like to apply a force to each atom along this normal vector.
Can I implement this kind of routine in a CHARMM code ?
Thank you in advance.


Edited by cola (08/10/05 01:31 PM)

Top
#7797 - 08/10/05 02:24 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: cola]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
It depends on what you mean by "CHARMM code". I doubt that forces relative to a re-orienting vector can be done via a CHARMM script using available commands. I'm certain it could be done in the Fortran source. Some care might be required to make sure other aspects of the code aren't broken, such as energy conservation.
_________________________
Rick Venable
computational chemist


Top
#7798 - 08/10/05 04:11 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: rmv]
cola Offline
Forum Member

Registered: 07/28/04
Posts: 113
Thank you for your reply.
I meant whether it can be done within a CHARMM script.
I want to keep the applying force normal to the plane always during the simulation.
Do you mean that it is not possible to find the positions ( and the normal vector) because there is no such kind of command in CHARMM ?

Top
#7799 - 08/10/05 04:56 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: cola]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
Correct, no existing command meets all of your requirements, esp. changing the direction of the applied force as the molecules rotates during dynamics.

The normal to a plane defined by a group of atom can be determined via

COOR LSQP NORM atom-selection

which will store the x,y,z components of the normal vector in the substitution values ?XAXI ?YAXI and ?ZAXI; see corman.doc

The PULL command (cons.doc) applies a force, but the vector orientation is fixed for the duration of the run.

Many of the commands in cons.doc are actually adding an applied force, but the vectors are typically based on atom coordinate positions in some way, as opposed to a derived vector (such as orthogonal to a plane formed by 3 atoms).
_________________________
Rick Venable
computational chemist


Top
#7800 - 08/10/05 07:01 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: rmv]
cola Offline
Forum Member

Registered: 07/28/04
Posts: 113
I'm a little confused.
If I know how to get the normal vector, then I think I can make a loop.

First, get a normal vector, and apply a force.
Then, after a few steps, get the normal vector again, and apply the force in a new direction.

Am I missing some points ?

Thank you in advance.

Top
#7801 - 08/10/05 09:28 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: cola]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
It could be done that way, but saving trajectories would be complicated, and the overhead of restarting dynamics would be fairly inefficient; it would be much, much slower than regular MD.

You could probably run some "normal" dynamics first, and get the P2 correlation function of a unit vector between two of the atoms that will define the plane. The time decay could be used to gauge a reasonable estimate of how many dynamics steps you can do in the loop before the vector orientation changes too much.
_________________________
Rick Venable
computational chemist


Top
#7802 - 08/11/05 01:57 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: rmv]
cola Offline
Forum Member

Registered: 07/28/04
Posts: 113
I see.
Thank you for your suggestions.

Top
#7803 - 08/17/05 02:25 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: rmv]
cola Offline
Forum Member

Registered: 07/28/04
Posts: 113
Please let me ask you one more question.

I assume that there are two normal vectors on the plane (going up and going down).
Which one is saved in the XAXI, YAXI, and ZAXI ?

I mean, If I have 1,2, and 3 atoms, and make a normal vector using COOR LSQM NORM (1,2,3 atoms), then the direction of normal vector is given by the clockwise or counter-clockwise rotation ?

Top
#7804 - 08/17/05 02:42 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: cola]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
It may be arbitrary, dependent on the numerical solution to the best fit plane; I'm not aware of any way to enforce a "right hand rule" based on the order of selected atoms or anything like that. You may have to test the vector and invert it if it points the wrong way.
_________________________
Rick Venable
computational chemist


Top
#7805 - 08/17/05 05:31 PM Re: Is this possible in CHARMM (adding a routine) ? [Re: rmv]
cola Offline
Forum Member

Registered: 07/28/04
Posts: 113
Thank you so much.
I'll check this.

Top
#7806 - 09/12/05 06:27 PM Angle between two least-square-planes with LSQP [Re: cola]
ilja Offline
Forum Member

Registered: 05/19/04
Posts: 34
Loc: La Jolla, California
Just thought that this might be useful for some who would want to compute angles between two mean planes, say for two stacking aromatic rings.

Because the plane equation as solved in CHARMM can have two solutions (defined up to inversion of the normal vector) straightforward computation of angles may not give consistent results.

This can be trivially addressed with the following script chunk that monitors the normal direction w.r.t. selected cross products generated by arbitrary three atoms from each of the two planes:


set pi = 3.14159265358979
define plane1 sele ... end
define plane2 sele ... end

define nc1 sele atom1_from_plane1 end
define nc2 sele atom2_from_plane1 end
define nc3 sele atom3_from_plane1 end

define pc1 sele atom1_from_plane2 end
define pc2 sele atom2_from_plane2 end
define pc3 sele atom3_from_plane3 end

open write card unit 22 name out_file.lsqp

!Your standard way of reading trajectories
traj query unit @istart
calc nframes = ?nfile * @pcs

traj iread @istart nunits @pcs nfile @nframes skip 1

set f 1

label tloop
traj read

format (f16.8)

coor axis sele nc1 end sele nc2 end
set vx1 = ?xaxi
set vy1 = ?yaxi
set vz1 = ?zaxi

coor axis sele nc1 end sele nc3 end
set vx2 = ?xaxi
set vy2 = ?yaxi
set vz2 = ?zaxi

calc cpx = @vy1*@vz2 - @vy2*@vz1
calc cpy = - @vx1*@vz2 + @vx2*@vz1
calc cpz = @vx1*@vy2 - @vx2*@vy1

coor lsqp norm mass sele plane1 end
set x1 = ?xaxi
set y1 = ?yaxi
set z1 = ?zaxi
set r1 = ?raxi
!now let's check if the vectors are in phase
calc dp = @cpx*@x1 + @cpy*@y1 + @cpz*@z1
if @dp .gt. 0 then goto skipinv1

set x1 = -@x1
set y1 = -@y1
set z1 = -@z1

label skipinv1

coor axis sele pc1 end sele pc2 end
set vx1 = ?xaxi
set vy1 = ?yaxi
set vz1 = ?zaxi

coor axis sele pc1 end sele pc3 end
set vx2 = ?xaxi
set vy2 = ?yaxi
set vz2 = ?zaxi

calc cpx = @vy1*@vz2 - @vy2*@vz1
calc cpy = - @vx1*@vz2 + @vx2*@vz1
calc cpz = @vx1*@vy2 - @vx2*@vy1

coor lsqp norm mass sele plane2 end
set x2 = ?xaxi
set y2 = ?yaxi
set z2 = ?zaxi
set r2 = ?raxi

!now let's check if the vectors are in phase
calc dp = @cpx*@x2 + @cpy*@y2 + @cpz*@z2
if @dp .gt. 0 then goto skipinv2

set x2 = -@x2
set y2 = -@y2
set z2 = -@z2

label skipinv2

calc angl = @x1*@x2 + @y1*@y2 + @z1*@z2
calc angl = @angl/@r1
calc angl = @angl/@r2
calc angl = acos(@angl)*180.0/@pi
calc angl2 = 180.0 - @angl

write title unit 22
* @f @angl @angl2
*

format

incr f by 1
if @f .le. @nframes then goto tloop

stop


Edited by ilja (09/13/05 12:13 PM)

Top
Page 1 of 2 1 2 >

Moderator:  chmgr, John Legato, petrella