[ESPResSo] Gay-Berne particles: Explanation of modifications

[Michael Winokur]

Hello,

So here are the details

1. Added parameter rinertia which allows for different moments of inertia about the principle axes. 

tcl syntax:   part $part_id rinertia $inertiax  $inertiay  $inertiaz

Requires #define ROTATIONAL_INERTIA */ in myconfig.h

Many other subroutines have been modified for proper interprocess communication between multiple nodes (untested)

2. Also added parameter quatu  (p->r.quatu[0] ,  p->r.quatu[1], p->r.quatu[2]) which is the (0,0,1) unit vector calculated from the quaternion.  This vector is needed in many places and so it seemed efficient to recalculate it just once after updating omega and the quaternion.  See rotation.c for the coding.  This parameter is also needs to be set in the interprocess communication area.

3. For the subtracted Gay-Berne 
tcl syntax:   incr bond_id; inter $bond_id subt_gb $reserved 4.70; set subt_CC  $bond_id
This system is isomorphic with the subt_lj but allows one to subtract it for connected near neighbor Gay-Berne pairs.
Note that the Gay-Berne syntax is unchanges
tcl syntax:  inter 1 1 gay-berne [expr 0.54*$ljmult] [expr (6.45)] 8. [expr 2.03/4.65] 3.67 2 1

4. For connected Gay-Berne particles:
tcl syntax:  inter $bond_id gb_dihedral 2 [expr -0.0923*0.5/$denom] [expr 2.*[PI]] 4 [expr -2.112*0.5/$denom] [expr 4.*[PI]] 6  [expr -0.3916*0.5/$denom] [expr 6.*[PI]] 0 -30.001000 0 0 -30.00100 0; set bond1_RCCR $bond_id

Here there are 15 parameters (5 sets of three).  The first three sets follow the updated dihedral interaction syntax and allow one to customize the torsional pair potential.  In my case I use ab-initio generated parameters published for biphenyl.  The forth set interaction specifies that the GB particle (0, 0 1) axis is perpendicular to the unit vector connecting the two GB particles.  In the case of biphenyl this keeps the plan of the benzene rings parallel to the bond connecting the two rings.  The fifth set interaction requires the particle submapping and fixes the azimuthal orientation of the GB particle.  In the case of biphenyl it keeps the two carbons making up the interring bond in from rotating in their respective planes.

5. Subparticle mapping
File "subparticle.xyz"
-----------------------------------------------------------------------
6
REMARK
C   -1.212   0.700   0.000
C   -1.212  -0.700   0.000
C    0.000  -1.400   0.000
C    1.212  -0.700   0.000
C    1.212   0.700   0.000
C   -0.000   1.400   0.000
2 3
0 1 2 2
0 2 3 3
------------------------------------------------------------------------- 
line 1  is the number of atoms in the subparticle.xyz file
line 2 is a comment line
lines 3-8  atom type and their x y z cooridinates
line 9  first integer is for the number mapping lines to read in and exact mapping needed for the fifth parameter set in gb_dihedral 
(i.e., which carbons form the C-C bond between two rings)
Note:  2 is for the number of subsequent terms, 3 is for the 0 to N X  0 to N array bounds of the map 
line 10 the first mapping between connected GB pairs
line 11 the second mapping between connected GB pairs
Four terms per line are necessary and, because the largest integer in any one line is 3 the second term in line 9 is 3.  
In the line "0 1 2 2" this signifies that the connected pair consists of GB particles 0 and 1 (or 1 and 0) using the carbon atom 2 x,y, z position (line 5) for GB particle 0 and the carbon atom position 2 for GB particle 1.  These allow for specifying unit vectors and don't reflect the two carbons forming the chemical bonds.  
In the line "0 2 3 3" this signifies that connected pair consists of  GB
particles 0 and 2 (or 2 and 0) using the carbon atom 23x,y, z position
(line 5) for GB particle 0 and the carbon atom position 3 for GB
particle 2.  Again, these allow for specifying unit vectors and don't reflect the two carbons forming the chemical bonds. 
In this case the vmd mapping and gb_inter interaction are using the same atoms but they could be independently specified.
Although a bit cumbersome (possibly confusing) this is the best I could do.  (Perhaps it would have been better to use another tcl parameter but that would mean I would have to further rewrite gb_inter.h to add all this feature).  Once imd.h and imd.c were modified, it was easier to use a separate data file.

The following is necessary for espresso to know about the subparticle mapping
It needs to be done early on (before configuring the vmd communications)
tcl syntax:        catch {imd subparticle "subparticle.xyz"} res''
The above tcl syntax is needed for the gb_dihed.c subroutine even if one does not plan to run the vmd software

Later in the tcl program this mapping needs to be further clarified
1 tcl syntax:  part $part_id sub_par 0 5 0
2 tcl syntax:  part $part_id sub_par 0 5 1
3 tcl syntax:  part $part_id sub_par 0 5 2
Here the first two parameters 0 & 5 specify the six carbon atoms given while the third if for the particular GB particle.  
 If there were 18 atoms in subparticle.xyz one could use  0 to 5 then 6 to 11 and finally 12 to 17. 
In the case of orthoterphenyl (i.e., three GB particles) the first term specifies a phenylene ring mapping (here carbons only) we reuse the same mapping three times for each of the three GB particles.  

This was not the easiest bit to program as, I have already noted, I am by no means a C-programming expert.

Please ask if you've any questions.

Cheers,

Michael