From f8a5e03a51c6aa29a2f9077e3387f102ef9344ec Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Fri, 20 Mar 2015 17:58:42 +0100
Subject: [PATCH] future releases of OpenMM will support triclinic boxes. We
have to make sure, that the angles get written properly into dcd files.
---
modules/mol/mm/src/mm_observer.cc | 29 +++++++++++++++++------------
1 file changed, 17 insertions(+), 12 deletions(-)
diff --git a/modules/mol/mm/src/mm_observer.cc b/modules/mol/mm/src/mm_observer.cc
index c8af30c12..aa0d5c5ee 100644
--- a/modules/mol/mm/src/mm_observer.cc
+++ b/modules/mol/mm/src/mm_observer.cc
@@ -130,18 +130,23 @@ void TrajWriter::Notify(){
// ucell
int32_t bsize=48; // ucell block size, 6 doubles
- //context_->getSystem().getDefaultPeriodicBoxVectors(ucell_a,ucell_b,ucell_c);
-
- // a,alpha,b,beta,gamma,c (don't ask)
- // OpenMM only supports rectangular boxes, since we have to add
- // the cosines of the angles, we can directly set the values to 0.0
-
- double ucell[]={10 * sqrt(ucell_a.dot(ucell_a)),
- 0.0,
- 10 * sqrt(ucell_b.dot(ucell_b)),
- 0.0,
- 0.0,
- 10 * sqrt(ucell_c.dot(ucell_c))};
+
+ double a_length = 10 * sqrt(ucell_a.dot(ucell_a)); //directly transform
+ double b_length = 10 * sqrt(ucell_b.dot(ucell_b)); //stuff to Angstrom
+ double c_length = 10 * sqrt(ucell_c.dot(ucell_c));
+
+ //note, that following angles have a compatibility issue
+ //CHARMM seems to write out the cosines of the actual
+ //unit cell angles, while NAMD seems to write the angles itself...
+ //we stick with the CHARMM way of doing things
+
+ double alpha = ucell_b.dot(ucell_c)/(b_length*c_length);
+ double beta = ucell_c.dot(ucell_a)/(c_length*a_length);
+ double gamma = ucell_a.dot(ucell_b)/(a_length*b_length);
+
+ double ucell[]={a_length, alpha, b_length,
+ beta, gamma, c_length};
+
stream_.write(reinterpret_cast<char*>(&bsize),4);
stream_.write(reinterpret_cast<char*>(ucell),bsize);
stream_.write(reinterpret_cast<char*>(&bsize),4);
--
GitLab