Commit ac481dd5 authored by Dr. Bogdan Tanygin's avatar Dr. Bogdan Tanygin Committed by GitHub

Merge pull request #145 from espressomd/python

Espresso head merge
parents 6ab6af86 272a3612
This diff is collapsed.
# Copyright (C) 2013,2014,2015,2016 The ESPResSo project
#
# This file is part of ESPResSo.
#
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
set lock ".lock-[pid]"
proc checklock {} {
global lock
if {! [file exists $lock]} {quit} {after 200 checklock}
}
menu main off
axes location off
mol load vtf case1.vtf
rock off
mol modstyle 0 0 CPK 1.500000 0.600000 8.000000 6.000000
mol modcolor 0 0 Name
exec touch $lock
while { [catch {imd connect localhost 10000}] } { after 500 }
checklock
This diff is collapsed.
This diff is collapsed.
# Copyright (C) 2013,2014,2015,2016 The ESPResSo project
#
# This file is part of ESPResSo.
#
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
set lock ".lock-[pid]"
proc checklock {} {
global lock
if {! [file exists $lock]} {quit} {after 200 checklock}
}
menu main off
axes location off
mol load vtf case2.vtf
rock off
rotate y by 90
scale by 2.1
mol modstyle 0 0 CPK 1.500000 0.600000 8.000000 6.000000
mol modcolor 0 0 Name
exec touch $lock
while { [catch {imd connect localhost 10000}] } { after 500 }
checklock
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
26 Axel (29.06.10 @ 17:39h)
41 Axel (29.06.10 @ 17:36h)
46 Axel (13.09.13 @ 09:40h)
47 Axel (12.09.13 @ 14:34h)
50 Axel (13.09.13 @ 09:37h)
60 Kurt BigBoss Kremer
This diff is collapsed.
This diff is collapsed.
......@@ -1213,8 +1213,9 @@ void force_and_velocity_display() {
/** @TODO: This needs to go!! */
int python_integrate(int n_steps, bool recalc_forces, bool reuse_forces) {
int python_integrate(int n_steps, bool recalc_forces, bool reuse_forces_par) {
int reuse_forces = 0;
reuse_forces = reuse_forces_par;
INTEG_TRACE(fprintf(stderr, "%d: integrate:\n", this_node));
if (recalc_forces) {
......
......@@ -152,9 +152,6 @@ cdef class Actor:
class Actors:
"""
.. todo:: implement def remove(self, actor)
"""
active_actors = []
......@@ -168,6 +165,15 @@ class Actors:
actor._activate()
else:
raise ThereCanOnlyBeOne(actor)
def remove(self, actor):
self._remove_actor(actor)
def _remove_actor(self, actor):
if not actor in self.active_actors:
raise Exception("Actor is not active")
actor._deactivate()
self.active_actors.remove(actor)
def __str__(self):
print("Active Actors:")
......@@ -187,7 +193,4 @@ class Actors:
def __delitem__(self, idx):
actor = self[idx]
if not actor in self.active_actors:
raise Exception("Actor is not active")
actor._deactivate()
self.active_actors.remove(actor)
self._remove_actor(actor)
......@@ -27,6 +27,11 @@ IF DIPOLES == 1:
int mdds_set_params(int n_cut)
int Ncut_off_magnetic_dipolar_direct_sum
IF (CUDA == 1) and (ROTATION == 1):
cdef extern from "actor/DipolarDirectSum.hpp":
void activate_dipolar_direct_sum_gpu()
void deactivate_dipolar_direct_sum_gpu()
IF DP3M == 1:
from p3m_common cimport p3m_parameter_struct
......
......@@ -63,6 +63,8 @@ IF DIPOLES == 1:
else:
self._params["bjerrum_length"] = self.params[
"prefactor"] / temperature
# also necessary on 1 CPU or GPU, does more than just broadcasting
mpi_bcast_coulomb_params()
def get_params(self):
self._params = self._get_params_from_es_core()
......@@ -72,6 +74,7 @@ IF DIPOLES == 1:
return coulomb.Dmethod
def _deactivate_method(self):
dipolar_set_Dbjerrum(0.0)
coulomb.Dmethod = DIPOLAR_NONE
mpi_bcast_coulomb_params()
......@@ -263,3 +266,33 @@ IF DIPOLES == 1:
if mdds_set_params(self._params["n_replica"]):
raise Exception(
"Could not activate magnetostatics method " + self.__class__.__name__)
IF (CUDA == 1) and (DIPOLES == 1) and (ROTATION == 1):
cdef class DipolarDirectSumGpu(MagnetostaticInteraction):
"""Calculates magnetostatic interactions by direct summation over all
pairs. If the system has periodic boundaries, the minimum image
convention is applied."""
def default_params(self):
return {}
def required_keys(self):
return ()
def valid_keys(self):
return ("bjerrum_length", "prefactor")
def _get_params_from_es_core(self):
return {"prefactor": coulomb.Dprefactor}
def _activate_method(self):
self._set_params_in_es_core()
def _deactivate_method(self):
super(type(self),self)._deactivate_method()
deactivate_dipolar_direct_sum_gpu()
def _set_params_in_es_core(self):
self.set_magnetostatics_prefactor()
activate_dipolar_direct_sum_gpu()
from espressomd.utils import to_char_pointer
from espressomd.utils import to_char_pointer, to_str
cdef class PScriptInterface:
def __init__(self, name=None, policy="GLOBAL"):
......@@ -29,7 +29,7 @@ cdef class PScriptInterface:
parameters = []
for p in self.sip.get().valid_parameters():
parameters.append(p.first)
parameters.append(to_str(p.first))
return parameters
......@@ -159,14 +159,14 @@ cdef class PScriptInterface:
raise Exception("Unkown type")
def get_parameter(self, name):
cdef Variant value = self.sip.get().get_parameter(name)
cdef Variant value = self.sip.get().get_parameter(to_char_pointer(name))
return self.variant_to_python_object(value)
def get_params(self):
cdef map[string, Variant] params = self.sip.get().get_parameters()
odict = {}
for pair in params:
odict[pair.first] = self.variant_to_python_object(pair.second)
odict[to_str(pair.first)] = self.variant_to_python_object(pair.second)
return odict
......@@ -184,7 +184,7 @@ class ScriptInterfaceHelper(PScriptInterface):
def __getattr__(self, attr):
if attr in self._valid_parameters():
return self.get_parameter(attr)
return self.get_parameter(to_char_pointer(attr))
else:
try:
return self.__dict__[attr]
......
......@@ -218,12 +218,19 @@ cdef class Thermostat:
"""
scalar_gamma_def = True
scalar_gamma_rot_def = True
IF PARTICLE_ANISOTROPY:
if isinstance(gamma, list):
scalar_gamma_def = False
else:
scalar_gamma_def = True
IF ROTATIONAL_INERTIA:
if isinstance(gamma_rotation, list):
scalar_gamma_rot_def = False
else:
scalar_gamma_rot_def = True
if kT is None or gamma is None:
raise ValueError(
"Both, kT and gamma have to be given as keyword args")
......@@ -234,6 +241,13 @@ cdef class Thermostat:
else:
utils.check_type_or_throw_except(
gamma, 3, float, "diagonal elements of the gamma tensor must be numbers")
if gamma_rotation is not None:
if scalar_gamma_rot_def:
utils.check_type_or_throw_except(
gamma_rotation, 1, float, "gamma_rotation must be a number")
else:
utils.check_type_or_throw_except(
gamma_rotation, 3, float, "diagonal elements of the gamma_rotation tensor must be numbers")
if scalar_gamma_def:
if float(kT) < 0. or float(gamma) < 0.:
......@@ -243,16 +257,16 @@ cdef class Thermostat:
if float(kT) < 0. or float(gamma[0]) < 0. or float(gamma[1]) < 0. or float(gamma[2]) < 0.:
raise ValueError(
"temperature and diagonal elements of the gamma tensor must be positive numbers")
global langevin_gamma_rotation
IF ROTATION:
if gamma_rotation is not None:
IF ROTATIONAL_INERTIA:
langevin_gamma_rotation[0] = gamma_rotation[0]
langevin_gamma_rotation[1] = gamma_rotation[1]
langevin_gamma_rotation[2] = gamma_rotation[2]
ELSE:
langevin_gamma_rotation = gamma_rotation
if gamma_rotation is not None:
if scalar_gamma_rot_def:
if float(gamma_rotation) < 0.:
raise ValueError(
"gamma_rotation must be positive number")
else:
if float(gamma_rotation[0]) < 0. or float(gamma_rotation[1]) < 0. or float(gamma_rotation[2]) < 0.:
raise ValueError(
"diagonal elements of the gamma_rotation tensor must be positive numbers")
global temperature
temperature = float(kT)
global langevin_gamma
......@@ -264,14 +278,52 @@ cdef class Thermostat:
else:
langevin_gamma[0] = gamma[0]
langevin_gamma[1] = gamma[1]
langevin_gamma[2] = gamma[3]
langevin_gamma[2] = gamma[2]
ELSE:
langevin_gamma = float(gamma)
global langevin_gamma_rotation
IF ROTATION:
if gamma_rotation is not None:
IF ROTATIONAL_INERTIA:
if scalar_gamma_rot_def:
langevin_gamma_rotation[0] = gamma_rotation
langevin_gamma_rotation[1] = gamma_rotation
langevin_gamma_rotation[2] = gamma_rotation
else:
langevin_gamma_rotation[0] = gamma_rotation[0]
langevin_gamma_rotation[1] = gamma_rotation[1]
langevin_gamma_rotation[2] = gamma_rotation[2]
ELSE:
if scalar_gamma_rot_def:
langevin_gamma_rotation = gamma_rotation
else:
raise ValueError(
"gamma_rotation must be a scalar since feature ROTATIONAL_INERTIA is disabled")
else:
IF ROTATIONAL_INERTIA:
IF PARTICLE_ANISOTROPY:
langevin_gamma_rotation[0] = langevin_gamma[0]
langevin_gamma_rotation[1] = langevin_gamma[1]
langevin_gamma_rotation[2] = langevin_gamma[2]
ELSE:
langevin_gamma_rotation[0] = langevin_gamma
langevin_gamma_rotation[1] = langevin_gamma
langevin_gamma_rotation[2] = langevin_gamma
ELSE:
IF PARTICLE_ANISOTROPY:
raise ValueError(
"gamma_rotation scalar parameter is required")
ELSE:
langevin_gamma_rotation = langevin_gamma
global thermo_switch
thermo_switch = (thermo_switch | THERMO_LANGEVIN)
mpi_bcast_parameter(FIELD_THERMO_SWITCH)
mpi_bcast_parameter(FIELD_TEMPERATURE)
mpi_bcast_parameter(FIELD_LANGEVIN_GAMMA)
IF ROTATION:
mpi_bcast_parameter(FIELD_LANGEVIN_GAMMA_ROTATION)
return True
IF LB_GPU or LB:
......
......@@ -3,6 +3,7 @@ set(py_tests bondedInteractions.py
constraint_shape_based.py
coulomb_cloud_wall.py
correlation.py
dawaanr-and-dds-gpu.py
electrostaticInteractions.py
engine_langevin.py
engine_lb.py
......
from __future__ import print_function
import unittest as ut
import numpy as np
import espressomd
from espressomd.interactions import *
from espressomd.magnetostatics import *
from espressomd.analyze import *
import math
from tests_common import *
from numpy import linalg as la
from numpy.random import random
from espressomd import has_features
@ut.skipIf(not has_features(["DIPOLES","CUDA","PARTIAL_PERIODIC","ROTATION"]),
"Features not available, skipping test!")
class DDSGPUTest(ut.TestCase):
longMessage = True
# Handle for espresso system
es = espressomd.System()
def vectorsTheSame(self,a,b):
tol = 3E-3
vec_len = la.norm(a - b)
if vec_len <= tol:
return True
else:
return False
def stopAll(self):
for i in range(len(self.es.part)):
self.es.part[i].v = np.array([0.0,0.0,0.0])
self.es.part[i].omega_body = np.array([0.0,0.0,0.0])
def run_test_case(self):
print("----------------------------------------------")
print("- Testcase dawaanr-and-dds-gpu.tcl")
print("----------------------------------------------")
pf_dds_gpu = 2.34
pf_dawaanr = 3.524
ratio_dawaanr_dds_gpu = pf_dawaanr / pf_dds_gpu
l = 15
self.es.box_l = [l, l, l]
self.es.periodicity = [0, 0, 0]
self.es.time_step = 1E-4
self.es.cell_system.skin = 0.1
part_dip = np.zeros((3))
for n in [ 110, 111, 540, 541 ]:
print("{0} particles".format(n))
dipole_modulus = 1.3
for i in range(n):
part_pos = np.array(random(3)) * l
costheta = 2 * random() - 1
sintheta = np.sin(np.arcsin(costheta))
phi = 2 * np.pi * random()
part_dip[0] = sintheta * np.cos(phi) * dipole_modulus
part_dip[1] = sintheta * np.sin(phi) * dipole_modulus
part_dip[2] = costheta * dipole_modulus
self.es.part.add(id = i, type = 0, pos = part_pos, dip = part_dip, v = np.array([0,0,0]), omega_body = np.array([0,0,0]))
self.es.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=10.0, sigma=0.5,
cutoff=0.55, shift="auto")
self.es.thermostat.set_langevin(kT=0.0, gamma=10.0)
self.es.integrator.set_steepest_descent(f_max=0.0, gamma=0.1, max_displacement=0.1)
self.es.integrator.run(500)
self.stopAll()
self.es.integrator.set_vv()
self.es.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=0.0, sigma=0.0,
cutoff=0.0, shift=0.0)
self.es.cell_system.skin = 0.0
self.es.time_step = 0.01
self.es.thermostat.turn_off()
# gamma should be zero in order to avoid the noise term in force and torque
self.es.thermostat.set_langevin(kT=1.297, gamma=0.0)
dds_cpu = DipolarDirectSumCpu(bjerrum_length = pf_dawaanr)
self.es.actors.add(dds_cpu)
self.es.integrator.run(steps = 0,recalc_forces = True)
dawaanr_f = []
dawaanr_t = []
for i in range(n):
dawaanr_f.append(self.es.part[i].f)
dawaanr_t.append(self.es.part[i].torque_lab)
dawaanr_e = Analysis(self.es).energy()["total"]
del dds_cpu
for i in range(len(self.es.actors.active_actors)):
self.es.actors.remove(self.es.actors.active_actors[i])
self.es.integrator.run(steps = 0,recalc_forces = True)
dds_gpu = DipolarDirectSumGpu(bjerrum_length = pf_dds_gpu)
self.es.actors.add(dds_gpu)
self.es.integrator.run(steps = 0,recalc_forces = True)
ddsgpu_f = []
ddsgpu_t = []
for i in range(n):
ddsgpu_f.append(self.es.part[i].f)
ddsgpu_t.append(self.es.part[i].torque_lab)
ddsgpu_e = Analysis(self.es).energy()["total"]
# compare
for i in range(n):
self.assertTrue(self.vectorsTheSame(np.array(dawaanr_t[i]),ratio_dawaanr_dds_gpu * np.array(ddsgpu_t[i])), \
msg = 'Torques on particle do not match. i={0} dawaanr_t={1} ratio_dawaanr_dds_gpu*ddsgpu_t={2}'.format(i,np.array(dawaanr_t[i]), ratio_dawaanr_dds_gpu * np.array(ddsgpu_t[i])))
self.assertTrue(self.vectorsTheSame(np.array(dawaanr_f[i]),ratio_dawaanr_dds_gpu * np.array(ddsgpu_f[i])), \
msg = 'Forces on particle do not match: i={0} dawaanr_f={1} ratio_dawaanr_dds_gpu*ddsgpu_f={2}'.format(i,np.array(dawaanr_f[i]), ratio_dawaanr_dds_gpu * np.array(ddsgpu_f[i])))
self.assertTrue(abs(dawaanr_e - ddsgpu_e * ratio_dawaanr_dds_gpu) <= 0.001, \
msg = 'Energies for dawaanr {0} and dds_gpu {1} do not match.'.format(dawaanr_e,ratio_dawaanr_dds_gpu * ddsgpu_e))
self.es.integrator.run(steps = 0,recalc_forces = True)
del dds_gpu
for i in range(len(self.es.actors.active_actors)):
self.es.actors.remove(self.es.actors.active_actors[i])
#for i in reversed(range(len(self.es.part))):
# self.es.part[i].delete()
self.es.part.clear()
def test(self):
if (self.es.cell_system.get_state()["n_nodes"] > 1):
print("NOTE: Ignoring testcase for n_nodes > 1")
else:
self.run_test_case()
if __name__ == '__main__':
print("Features: ", espressomd.features())
ut.main()
......@@ -6,7 +6,7 @@ import math
@ut.skipIf(not espressomd.has_features(["MASS","ROTATIONAL_INERTIA"]),
"Features not available, skipping test!")
class ThermoTest(ut.TestCase):
class RotationalInertia(ut.TestCase):
longMessage = True
# Handle for espresso system
es = espressomd.System()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment