From b2b465571e9187c95f28348e9634809e82cb981a Mon Sep 17 00:00:00 2001
From: pinyili2 <pinyili2@illinois.edu>
Date: Tue, 8 Oct 2024 15:34:47 -0500
Subject: [PATCH] make mrdna competitable with current arbdmodel
---
mrdna/RELEASE-VERSION | 2 +-
mrdna/__init__.py | 2 +-
mrdna/marbdmodel/coords.py | 3 ++-
mrdna/marbdmodel/coords_arbd.py | 1 +
mrdna/segmentmodel.py | 18 ++++++++++++++----
5 files changed, 19 insertions(+), 7 deletions(-)
diff --git a/mrdna/RELEASE-VERSION b/mrdna/RELEASE-VERSION
index 12166bd..59d9bd8 100644
--- a/mrdna/RELEASE-VERSION
+++ b/mrdna/RELEASE-VERSION
@@ -1 +1 @@
-1.0a.dev171
+1.0a.dev177
diff --git a/mrdna/__init__.py b/mrdna/__init__.py
index a60c9d3..7eceda7 100644
--- a/mrdna/__init__.py
+++ b/mrdna/__init__.py
@@ -35,7 +35,7 @@ def get_resource_path(relative_path):
try:
from arbdmodel.coords import read_arbd_coordinates
except:
- from marbdmodel.coords import readArbdCoords
+ from .marbdmodel.coords import readArbdCoords
from .segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
from .simulate import multiresolution_simulation
from .model.dna_sequence import read_sequence_file
diff --git a/mrdna/marbdmodel/coords.py b/mrdna/marbdmodel/coords.py
index df1f145..d9837b0 100644
--- a/mrdna/marbdmodel/coords.py
+++ b/mrdna/marbdmodel/coords.py
@@ -110,6 +110,7 @@ def quaternion_to_matrix(q):
return np.array(R)
def quaternion_from_matrix( R ):
+ R=R.T
q = np.empty(4)
if R[2,2] < 0:
if R[0,0] > R[1,1]:
@@ -204,7 +205,7 @@ def quaternion_product(a, b):
def quaternion_inverse(q):
assert(len(q) == 4)
-quaternion_from_matrix( R ) qinv = np.array(q)
+ qinv = np.array(q)
qinv[1:] = -qinv[1:]
qinv = qinv / (q[0]**2 + q[1:].dot(q[1:]))
return qinv
diff --git a/mrdna/marbdmodel/coords_arbd.py b/mrdna/marbdmodel/coords_arbd.py
index c9ecd7f..2df61ed 100644
--- a/mrdna/marbdmodel/coords_arbd.py
+++ b/mrdna/marbdmodel/coords_arbd.py
@@ -111,6 +111,7 @@ def quaternion_to_matrix(q):
return np.array(R)
def quaternion_from_matrix( R ):
+ R=R.T
q = np.empty(4)
if R[2,2] < 0:
if R[0,0] > R[1,1]:
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index d09c45b..e77f8d0 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -524,7 +524,7 @@ class Segment(ConnectableElement, Group):
# orientation_bond = HarmonicBond(10,2)
if sys_arbd==False:
- orientation_bond = HarmonicBond(30,1.5, range_ = (0,500) )
+ orientation_bond = HarmonicBond(30,1.5, rRange = (0,500) )
else:
orientation_bond = HarmonicBond(30,1.5, range_ = (0,500) )
ssDNA_particle = ParticleType("S",
@@ -2028,9 +2028,16 @@ class SegmentModel(ArbdModel):
if key not in self._bonded_potential:
assert( kSpring >= 0 )
if type_ == "bond":
- self._bonded_potential[key] = HarmonicBond(kSpring,d, range_=(0,1200), max_potential=max_potential)
+ if sys_arbd==True:
+ self._bonded_potential[key] = HarmonicBond(kSpring,d, range_=(0,1200), max_potential=max_potential)
+ else:
+ self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,1200), max_potential=max_potential)
elif type_ == "gbond":
- self._bonded_potential[key] = HarmonicBond(kSpring,d, range_=(0,1200), max_potential=max_potential, correct_geometry=True, temperature=self.temperature)
+ if sys_arbd==True:
+ self._bonded_potential[key] = HarmonicBond(kSpring,d, range_=(0,1200), max_potential=max_potential, correct_geometry=True, temperature=self.temperature)
+ else:
+ self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,1200), max_potential=max_potential, correct_geometry=True, temperature=self.temperature)
+
elif type_ == "angle":
self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
# , resolution = 1, maxForce=0.1)
@@ -2061,7 +2068,10 @@ class SegmentModel(ArbdModel):
assert( d > 0.2 )
if key not in self._bonded_potential:
kT = self.temperature * 0.0019872065 # kcal/mol
- self._bonded_potential[key] = WLCSKBond( d, lp, kT, range_=(0,1200) )
+ if sys_arbd is True:
+ self._bonded_potential[key] = WLCSKBond( d, lp, kT, range_=(0,1200) )
+ else:
+ self._bonded_potential[key] = WLCSKBond( d, lp, kT, rRange=(0,1200) )
return self._bonded_potential[key]
def _get_wlc_sk_angle_potential(self, d):
--
GitLab