Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
Is this possible in CHARMM (adding a routine) ?
#7796 08/10/05 04:12 PM
Joined: Jul 2004
Posts: 113
C
cola Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Jul 2004
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.

Last edited by cola; 08/10/05 05:31 PM.
Re: Is this possible in CHARMM (adding a routine) ?
cola #7797 08/10/05 06:24 PM
Joined: Sep 2003
Posts: 8,450
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,450
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

Re: Is this possible in CHARMM (adding a routine) ?
rmv #7798 08/10/05 08:11 PM
Joined: Jul 2004
Posts: 113
C
cola Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Jul 2004
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 ?

Re: Is this possible in CHARMM (adding a routine) ?
cola #7799 08/10/05 08:56 PM
Joined: Sep 2003
Posts: 8,450
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,450
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

Re: Is this possible in CHARMM (adding a routine) ?
rmv #7800 08/10/05 11:01 PM
Joined: Jul 2004
Posts: 113
C
cola Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Jul 2004
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.

Re: Is this possible in CHARMM (adding a routine) ?
cola #7801 08/11/05 01:28 AM
Joined: Sep 2003
Posts: 8,450
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,450
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

Re: Is this possible in CHARMM (adding a routine) ?
rmv #7802 08/11/05 05:57 PM
Joined: Jul 2004
Posts: 113
C
cola Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Jul 2004
Posts: 113
I see.
Thank you for your suggestions.

Re: Is this possible in CHARMM (adding a routine) ?
rmv #7803 08/17/05 06:25 PM
Joined: Jul 2004
Posts: 113
C
cola Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Jul 2004
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 ?

Re: Is this possible in CHARMM (adding a routine) ?
cola #7804 08/17/05 06:42 PM
Joined: Sep 2003
Posts: 8,450
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,450
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

Re: Is this possible in CHARMM (adding a routine) ?
rmv #7805 08/17/05 09:31 PM
Joined: Jul 2004
Posts: 113
C
cola Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Jul 2004
Posts: 113
Thank you so much.
I'll check this.

Angle between two least-square-planes with LSQP
cola #7806 09/12/05 10:27 PM
Joined: May 2004
Posts: 34
I
Forum Member
Offline
Forum Member
I
Joined: May 2004
Posts: 34
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

Last edited by ilja; 09/13/05 04:13 PM.
Page 1 of 2 1 2

Moderated by  chmgr, John Legato, petrella 

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.009s Queries: 36 (0.004s) Memory: 0.9965 MB (Peak: 1.1388 MB) Data Comp: Off Server Time: 2020-06-03 20:32:18 UTC
Valid HTML 5 and Valid CSS