Sorry for the long delay in communications but I encountered many unforeseen problems in suitably modifying the code. Basically the original paper I was using for the equations of motion was broken (which I failed to initially recognize) and there were limitations in the original espresso GB source (I'm not convinced it was 100% correct). It wasn't until I looked at another paper and the source code from IMD (The ITAP Molecular Dynamics Program) that things started to work properly. [In the later I believe there is a small but significant sign error in their dU/dc (i.e., c is the torsion angle between GB particles) term.] I'm not a theorist and so I am not comfortable with tackling these kind of partial derivative.
Here is what I have added:
1. Gay-Berne particle pair interaction gb.c Espresso typically allows for an r cutoff which I've added to my version of gb.c but the IMD source doesn't include this option.
2. Gay-Berne subtracted interaction subt_gb.c 3. connected GB particles gb_dihed.c 4. Anisotropic Rotational inertia parameters. 5. A submapping scheme so that vmd program can display a subparticle structure.
other subroutines modified in the process energy.h, energy.c, forces.h, forces.c, communication.h, communication.c, thermostat.c, thermostat.h (I'm not 100% certain that I really reproduce the results dictated by the equipartition theorem.) particle_data.c, particle_data.h, interaction_data.c, interaction_data.h, imd.c, imd.h, rotation.c and myconfig.h
At the moment the GB interaction is between identical particles. It should be relatively easy to introduce a scaling factor that allows for different shape objects (e.g. the interaction between a GB and united atom pair ) but I need to study the MD of ortho-terphenyl first before I start venturing into this area.
It is really is nice to see F=ma and T=I alpha...and at zero temperature the runs downhill.