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