From 1e1f9ac232cfaef3cca71cd1a52506969c52e2b4 Mon Sep 17 00:00:00 2001
From: Carl Poelking <cp605@cam.ac.uk>
Date: Thu, 21 Sep 2017 12:39:49 +0100
Subject: [PATCH] Harmonic net, TODO remove com motion.
---
src/soap/soapy/dimred.py | 13 +++++++++++--
test/test_harmonicnet/network.py | 23 +++++++++++++++++++++++
2 files changed, 34 insertions(+), 2 deletions(-)
create mode 100644 test/test_harmonicnet/network.py
diff --git a/src/soap/soapy/dimred.py b/src/soap/soapy/dimred.py
index 17e16e9..078de10 100644
--- a/src/soap/soapy/dimred.py
+++ b/src/soap/soapy/dimred.py
@@ -72,7 +72,7 @@ class BondNetwork(object):
beta = 1.
k_t = 0.85 # 0.85 # 0.85
self.H = beta*0.5*(
- 1.+np.tanh(alpha*(K-k_t))
+ 1.+np.tanh(alpha*(self.K-k_t))
)
# Universal force constant?
self.H.fill(0.28) # NOTE HACK
@@ -169,9 +169,9 @@ class BondNetwork(object):
- self.R0[self.idcs_noexcl][:,self.idcs_noexcl])**2)
else:
E = np.sum(0.5*self.H*(dR - self.R0)**2)
- print "Step %d %+1.7e/%+1.7e %+1.7e" % (i, rms, rms_max, E)
# Trajectory
if i % dn_out == 0:
+ print "Step %5d rms/max %+1.7e/%+1.7e energy %+1.7e" % (i, rms, rms_max, E)
scale = 100.
ofs.write('%d\n\n' % (self.N))
for i in range(self.N):
@@ -183,8 +183,17 @@ class BondNetwork(object):
if rms < rms_cut and rms_max < rms_cut:
print "*** converged ***"
break
+ print "Step %5d rms/max %+1.7e/%+1.7e energy %+1.7e" % (i, rms, rms_max, E)
ofs.close()
return self.X, self.Y
+ def write_confout(self, filename):
+ ofs = open(filename, 'w')
+ for i in range(self.N):
+ ofs.write('%10s %1.2e %+1.4e %+1.4e\n' % (
+ self.tags[i], self.M[i], self.X[i], self.Y[i]))
+ ofs.close()
+ return
+
def dimred_matrix(method, kmat=None, distmat=None, outfile=None, ix=None, symmetrize=False, prj_dimension=2):
if symmetrize:
diff --git a/test/test_harmonicnet/network.py b/test/test_harmonicnet/network.py
new file mode 100644
index 0000000..91992b6
--- /dev/null
+++ b/test/test_harmonicnet/network.py
@@ -0,0 +1,23 @@
+import numpy as np
+import soap
+import json
+
+# Nodes
+N = 100
+tags = [ 'P%03d' for i in range(N) ]
+kernel = 0.5*(1. + np.random.uniform(size=N*N).reshape((N,N)))
+masses = 10*np.random.uniform(size=N) + 1.
+
+# Network
+network = soap.soapy.BondNetwork(tags, kernel, masses)
+network.initialise(method='kernelpca')
+network.integrate_md(
+ n_steps=5000,
+ rms_cut=1e-7,
+ dt=0.01,
+ dn_out=100,
+ append_traj=False)
+
+# Write final frame
+network.write_confout('confout.txt')
+
--
GitLab