1# -*- coding: utf-8 -*- 2# MolMod is a collection of molecular modelling tools for python. 3# Copyright (C) 2007 - 2019 Toon Verstraelen <Toon.Verstraelen@UGent.be>, Center 4# for Molecular Modeling (CMM), Ghent University, Ghent, Belgium; all rights 5# reserved unless otherwise stated. 6# 7# This file is part of MolMod. 8# 9# MolMod is free software; you can redistribute it and/or 10# modify it under the terms of the GNU General Public License 11# as published by the Free Software Foundation; either version 3 12# of the License, or (at your option) any later version. 13# 14# MolMod is distributed in the hope that it will be useful, 15# but WITHOUT ANY WARRANTY; without even the implied warranty of 16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17# GNU General Public License for more details. 18# 19# You should have received a copy of the GNU General Public License 20# along with this program; if not, see <http://www.gnu.org/licenses/> 21# 22# -- 23 24 25from __future__ import division 26 27from builtins import range 28import numpy as np 29import pkg_resources 30 31import molmod.ic as ic 32from molmod import * 33 34 35num = 10 36big = 5 37small = 1 38eps = 1e-5 39 40 41def iter_bonds(): 42 counter = 0 43 while counter < num: 44 p1 = np.random.normal(0,big,3) 45 p2 = np.random.normal(0,big,3) 46 if np.linalg.norm(p1-p2) < small: 47 continue 48 yield p1, p2 49 counter += 1 50 51 52def iter_bends(): 53 counter = 0 54 while counter < num: 55 p1 = np.random.normal(0,big,3) 56 p2 = np.random.normal(0,big,3) 57 if np.linalg.norm(p1-p2) < small: 58 continue 59 p3 = np.random.normal(0,big,3) 60 if np.linalg.norm(p1-p3) < small: 61 continue 62 if np.linalg.norm(p2-p3) < small: 63 continue 64 yield p1, p2, p3 65 counter += 1 66 67 68def iter_diheds(): 69 counter = 0 70 while counter < num: 71 p1 = np.random.normal(0,big,3) 72 p2 = np.random.normal(0,big,3) 73 if np.linalg.norm(p1-p2) < small: 74 continue 75 d = (p1-p2)/np.linalg.norm(p1-p2) 76 p3 = np.random.normal(0,big,3) 77 p3_ortho = p3 - d*np.dot(p3-p1,d) 78 if np.linalg.norm(p3_ortho) < small: 79 continue 80 p4 = np.random.normal(0,big,3) 81 p4_ortho = p4 - d*np.dot(p4-p2,d) 82 if np.linalg.norm(p4_ortho) < small: 83 continue 84 yield p1, p2, p3, p4 85 counter += 1 86 87 88def iter_pairs(): 89 for counter in range(num): 90 while True: 91 x = np.random.normal(0,big,2) 92 if x[-1] > small: break 93 yield x 94 95 96def check_diff_ic(icfn, iterp, shape=(-1,3), period=None): 97 def fnv(x0, do_gradient=False): 98 q, g = icfn(x0.reshape(shape),1) 99 if do_gradient: 100 return q, g.ravel() 101 else: 102 return q 103 104 def fng(x0, do_gradient=False): 105 q, g, h = icfn(x0.reshape(shape),2) 106 if do_gradient: 107 return g.ravel(), h.reshape(g.size,g.size) 108 else: 109 return g.ravel() 110 111 for ps in iterp(): 112 x0 = np.array(ps).ravel() 113 dxs = np.random.normal(0, eps, (100, len(x0))) 114 check_delta(fnv, x0, dxs, period) 115 check_delta(fng, x0, dxs, period) 116 117 118def test_diff_bond(): 119 check_diff_ic(ic.bond_length, iter_bonds) 120 121 122def test_diff_dot(): 123 def my_dot(rs,deriv): 124 x = ic.Vector3(6,deriv,rs[0],(0,1,2)) 125 y = ic.Vector3(6,deriv,rs[1],(3,4,5)) 126 return ic.dot(x,y).results() 127 check_diff_ic(my_dot, iter_bonds) 128 129 130class MyCross(object): 131 def __init__(self, index): 132 self.index = index 133 134 def __call__(self, rs, deriv): 135 x = ic.Vector3(6,deriv,rs[0],(0,1,2)) 136 y = ic.Vector3(6,deriv,rs[1],(3,4,5)) 137 if self.index == 0: 138 return ic.cross(x,y).x.results() 139 elif self.index == 1: 140 return ic.cross(x,y).y.results() 141 elif self.index == 2: 142 return ic.cross(x,y).z.results() 143 else: 144 raise NotImplementedError 145 146 147def test_diff_cross(): 148 check_diff_ic(MyCross(0), iter_bonds) 149 check_diff_ic(MyCross(1), iter_bonds) 150 check_diff_ic(MyCross(2), iter_bonds) 151 152 153def test_diff_idiv(): 154 def my_div(t,deriv): 155 t0 = ic.Scalar(2,deriv,t[0],0) 156 t1 = ic.Scalar(2,deriv,t[1],1) 157 t0 /= t1 158 return t0.results() 159 check_diff_ic(my_div, iter_pairs, (2,)) 160 161 162def test_diff_bend_angle(): 163 check_diff_ic(ic.bend_angle, iter_bends) 164 165 166def test_diff_bend_cos(): 167 check_diff_ic(ic.bend_cos, iter_bends) 168 169 170def test_diff_imul(): 171 def my_mul(t,deriv): 172 t0 = ic.Scalar(2,deriv,t[0],0) 173 t1 = ic.Scalar(2,deriv,t[1],1) 174 t0 *= t1 175 return t0.results() 176 check_diff_ic(my_mul, iter_pairs, (2,)) 177 178 179def test_diff_dihed_cos(): 180 check_diff_ic(ic.dihed_cos, iter_diheds) 181 182 183def test_diff_dihed_angle(): 184 check_diff_ic(ic.dihed_angle, iter_diheds, period=2*np.pi) 185 186 187def iter_diheds_special(): 188 yield np.array([[+1,0,1], [0,0,0], [0,0,1], [1,0,1]]) 189 yield np.array([[0,+1,1], [0,0,0], [0,0,1], [1,0,1]]) 190 yield np.array([[-1,0,1], [0,0,0], [0,0,1], [1,0,1]]) 191 yield np.array([[0,-1,1], [0,0,0], [0,0,1], [1,0,1]]) 192 193 194def test_diff_dihed_angle_special(): 195 check_diff_ic(ic.dihed_angle, iter_diheds_special, period=2*np.pi) 196 197 198def test_diff_opbend_cos(): 199 check_diff_ic(ic.opbend_cos, iter_diheds) 200 201 202def test_diff_opbend_angle(): 203 check_diff_ic(ic.opbend_angle, iter_diheds) 204 205def test_diff_opbend_mcos(): 206 check_diff_ic(ic.opbend_mcos, iter_diheds) 207 208def test_diff_opbend_mangle(): 209 check_diff_ic(ic.opbend_mangle, iter_diheds) 210 211def test_diff_opbend_dist(): 212 check_diff_ic(ic.opbend_dist, iter_diheds) 213 214 215def test_dihedral_ethene(): 216 mol = Molecule.from_file(pkg_resources.resource_filename("molmod", "data/test/ethene.xyz")) 217 c = mol.coordinates.copy() 218 assert abs(ic.dihed_cos([c[2], c[0], c[3], c[5]])[0] - 1.0) < 1e-5 219 assert abs(ic.dihed_angle([c[2], c[0], c[3], c[5]])[0]) < 1e-5 220 for i in range(1000): 221 angle = np.random.uniform(-np.pi, np.pi) 222 radius = np.random.uniform(0, 5*angstrom) 223 offset = np.random.uniform(0, 5*angstrom) 224 c[5] = [ 225 offset, 226 -radius*np.cos(angle), 227 -radius*np.sin(angle), 228 ] 229 assert abs(ic.dihed_cos([c[2], c[0], c[3], c[5]])[0] - np.cos(angle)) < 1e-5 230 assert abs(ic.dihed_angle([c[2], c[0], c[3], c[5]])[0] - angle) < 1e-5 231 232 233def test_opbend_ethene(): 234 mol = Molecule.from_file(pkg_resources.resource_filename("molmod", "data/test/ethene.xyz")) 235 c = mol.coordinates.copy() 236 assert abs(ic.opbend_cos([c[0], c[5], c[4], c[3]])[0] - 1.0) < 1e-5 237 assert abs(ic.opbend_angle([c[0], c[5], c[4], c[3]])[0]) < 1e-5 238 assert abs(ic.opbend_mcos([c[0], c[5], c[4], c[3]])[0] - 1.0) < 1e-5 239 assert abs(ic.opbend_mangle([c[0], c[5], c[4], c[3]])[0]) < 1e-5 240 assert abs(ic.opbend_dist([c[0], c[5], c[4], c[3]])[0]) < 1e-5 241 for i in range(1000): 242 angle = np.random.uniform(-np.pi/2, np.pi/2) 243 radius = np.random.uniform(0, 5*angstrom) 244 #offset = np.random.uniform(0, 5*angstrom) 245 c[3] = [ 246 radius*np.cos(angle), 247 0.0, 248 radius*np.sin(angle), 249 ] 250 assert abs(ic.opbend_cos([c[0], c[5], c[4], c[3]])[0] - np.cos(angle)) < 1e-5 251 assert abs(ic.opbend_angle([c[0], c[5], c[4], c[3]])[0] - angle) < 1e-5 252