Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
Scientific Software
mmtk
Commits
3e132477
Commit
3e132477
authored
Dec 11, 2017
by
eric pellegrini
Browse files
removed oldnumeric deps from ChargeFit, Collections, ConfigIO and DCD
parent
3e974dbd
Changes
4
Hide whitespace changes
Inline
Side-by-side
MMTK/ChargeFit.py
View file @
3e132477
...
...
@@ -26,9 +26,10 @@ See also Examples/Miscellaneous/charge_fit.py.
__docformat__
=
'restructuredtext'
import
numpy
from
MMTK
import
Random
,
Units
,
Utility
from
Scientific.Geometry
import
Vector
from
Scientific
import
N
from
Scientific
import
LA
...
...
@@ -68,8 +69,8 @@ class ChargeFit(object):
if
npoints
<
natoms
:
raise
ValueError
(
"Not enough data points for fit"
)
m
=
N
.
zeros
((
npoints
,
natoms
),
N
.
F
loat
)
phi
=
N
.
zeros
((
npoints
,),
N
.
F
loat
)
m
=
numpy
.
zeros
((
npoints
,
natoms
),
numpy
.
f
loat
)
phi
=
numpy
.
zeros
((
npoints
,),
numpy
.
f
loat
)
i
=
0
for
conf
,
pointlist
in
points
.
items
():
for
r
,
p
in
pointlist
:
...
...
@@ -83,22 +84,22 @@ class ChargeFit(object):
phi_test
=
phi
if
constraints
is
not
None
:
phi
-=
N
.
dot
(
m
,
constraints
.
bi_c
)
m
=
N
.
dot
(
m
,
constraints
.
p
)
phi
-=
numpy
.
dot
(
m
,
constraints
.
bi_c
)
m
=
numpy
.
dot
(
m
,
constraints
.
p
)
c_rank
=
constraints
.
rank
else
:
c_rank
=
0
u
,
s
,
vt
=
LA
.
singular_value_decomposition
(
m
)
s_test
=
s
[:
len
(
s
)
-
c_rank
]
cutoff
=
1.e-10
*
N
.
maximum
.
reduce
(
s_test
)
nonzero
=
N
.
repeat
(
s_test
,
N
.
not_equal
(
s_test
,
0.
))
cutoff
=
1.e-10
*
numpy
.
maximum
.
reduce
(
s_test
)
nonzero
=
numpy
.
repeat
(
s_test
,
numpy
.
not_equal
(
s_test
,
0.
))
self
.
rank
=
len
(
nonzero
)
self
.
condition
=
N
.
maximum
.
reduce
(
nonzero
)
/
\
N
.
minimum
.
reduce
(
nonzero
)
self
.
effective_rank
=
N
.
add
.
reduce
(
N
.
greater
(
s
,
cutoff
))
self
.
condition
=
numpy
.
maximum
.
reduce
(
nonzero
)
/
\
numpy
.
minimum
.
reduce
(
nonzero
)
self
.
effective_rank
=
numpy
.
add
.
reduce
(
numpy
.
greater
(
s
,
cutoff
))
if
self
.
effective_rank
<
self
.
rank
:
self
.
effective_condition
=
N
.
maximum
.
reduce
(
nonzero
)
/
cutoff
self
.
effective_condition
=
numpy
.
maximum
.
reduce
(
nonzero
)
/
cutoff
else
:
self
.
effective_condition
=
self
.
condition
if
self
.
effective_rank
<
natoms
-
c_rank
:
...
...
@@ -110,15 +111,15 @@ class ChargeFit(object):
s
[
i
]
=
1.
/
s
[
i
]
else
:
s
[
i
]
=
0.
q
=
N
.
dot
(
N
.
transpose
(
vt
),
s
*
N
.
dot
(
N
.
transpose
(
u
)[:
natoms
,
:],
phi
))
q
=
numpy
.
dot
(
numpy
.
transpose
(
vt
),
s
*
numpy
.
dot
(
numpy
.
transpose
(
u
)[:
natoms
,
:],
phi
))
if
constraints
is
not
None
:
q
=
constraints
.
bi_c
+
N
.
dot
(
constraints
.
p
,
q
)
q
=
constraints
.
bi_c
+
numpy
.
dot
(
constraints
.
p
,
q
)
deviation
=
N
.
dot
(
m_test
,
q
)
-
phi_test
self
.
rms_error
=
N
.
sqrt
(
N
.
dot
(
deviation
,
deviation
))
deviation
=
N
.
fabs
(
deviation
/
phi_test
)
self
.
relative_rms_error
=
N
.
sqrt
(
N
.
dot
(
deviation
,
deviation
))
deviation
=
numpy
.
dot
(
m_test
,
q
)
-
phi_test
self
.
rms_error
=
numpy
.
sqrt
(
numpy
.
dot
(
deviation
,
deviation
))
deviation
=
numpy
.
fabs
(
deviation
/
phi_test
)
self
.
relative_rms_error
=
numpy
.
sqrt
(
numpy
.
dot
(
deviation
,
deviation
))
self
.
charges
=
{}
for
i
in
range
(
natoms
):
...
...
@@ -192,8 +193,8 @@ class ChargeConstraintSet(object):
self
.
atoms
=
atoms
natoms
=
len
(
self
.
atoms
)
nconst
=
sum
([
len
(
c
)
for
c
in
constraints
])
b
=
N
.
zeros
((
nconst
,
natoms
),
N
.
F
loat
)
c
=
N
.
zeros
((
nconst
,),
N
.
F
loat
)
b
=
numpy
.
zeros
((
nconst
,
natoms
),
numpy
.
f
loat
)
c
=
numpy
.
zeros
((
nconst
,),
numpy
.
f
loat
)
i
=
0
for
cons
in
constraints
:
cons
.
setCoefficients
(
self
.
atoms
,
b
,
c
,
i
)
...
...
@@ -205,11 +206,11 @@ class ChargeConstraintSet(object):
self
.
rank
=
self
.
rank
+
1
self
.
b
=
b
self
.
bi
=
LA
.
generalized_inverse
(
b
)
self
.
p
=
N
.
identity
(
natoms
)
-
N
.
dot
(
self
.
bi
,
self
.
b
)
self
.
p
=
numpy
.
identity
(
natoms
)
-
numpy
.
dot
(
self
.
bi
,
self
.
b
)
self
.
c
=
c
self
.
bi_c
=
N
.
dot
(
self
.
bi
,
c
)
c_test
=
N
.
dot
(
self
.
b
,
self
.
bi_c
)
if
N
.
add
.
reduce
((
c_test
-
c
)
**
2
)
/
nconst
>
1.e-12
:
self
.
bi_c
=
numpy
.
dot
(
self
.
bi
,
c
)
c_test
=
numpy
.
dot
(
self
.
b
,
self
.
bi_c
)
if
numpy
.
add
.
reduce
((
c_test
-
c
)
**
2
)
/
nconst
>
1.e-12
:
Utility
.
warning
(
"The charge constraints are inconsistent."
" They will be applied as a least-squares"
" condition."
)
...
...
MMTK/Collections.py
View file @
3e132477
...
...
@@ -9,11 +9,12 @@ Collections of chemical objects
__docformat__
=
'restructuredtext'
import
numpy
from
MMTK
import
Utility
,
Units
,
ParticleProperties
,
Visualization
from
MMTK.Geometry
import
superpositionFit
from
Scientific.Geometry
import
Vector
,
Tensor
,
Objects3D
from
Scientific.Geometry
import
Transformation
from
Scientific
import
N
import
copy
,
itertools
,
types
#
...
...
@@ -161,7 +162,7 @@ class GroupOfAtoms(object):
"""
com
,
i
=
self
.
centerAndMomentOfInertia
(
conf
)
pmi
=
i
.
eigenvalues
()
return
N
.
sort
(
Units
.
h
/
(
8.
*
N
.
pi
*
N
.
pi
*
pmi
))[::
-
1
]
return
numpy
.
sort
(
Units
.
h
/
(
8.
*
numpy
.
pi
*
numpy
.
pi
*
pmi
))[::
-
1
]
def
boundingBox
(
self
,
conf
=
None
):
"""
...
...
@@ -178,8 +179,8 @@ class GroupOfAtoms(object):
max
=
min
for
a
in
atoms
[
1
:]:
r
=
a
.
position
(
conf
).
array
min
=
N
.
minimum
(
min
,
r
)
max
=
N
.
maximum
(
max
,
r
)
min
=
numpy
.
minimum
(
min
,
r
)
max
=
numpy
.
maximum
(
max
,
r
)
return
Vector
(
min
),
Vector
(
max
)
def
boundingSphere
(
self
,
conf
=
None
):
...
...
@@ -220,7 +221,7 @@ class GroupOfAtoms(object):
dr
=
universe
.
distanceVector
(
a
.
position
(
conf1
),
a
.
position
(
conf2
))
m
+=
ma
rms
+=
ma
*
dr
*
dr
return
N
.
sqrt
(
rms
/
m
)
return
numpy
.
sqrt
(
rms
/
m
)
def
findTransformationAsQuaternion
(
self
,
conf1
,
conf2
=
None
):
universe
=
self
.
universe
()
...
...
@@ -325,20 +326,20 @@ class GroupOfAtoms(object):
if
determinant
(
diag
.
array
)
<
0
:
diag
.
array
[
0
]
=
-
diag
.
array
[
0
]
if
repr
!=
None
:
seq
=
N
.
argsort
(
ev
)
seq
=
numpy
.
argsort
(
ev
)
if
repr
==
'Ir'
:
seq
=
N
.
array
([
seq
[
1
],
seq
[
2
],
seq
[
0
]])
seq
=
numpy
.
array
([
seq
[
1
],
seq
[
2
],
seq
[
0
]])
elif
repr
==
'IIr'
:
seq
=
N
.
array
([
seq
[
2
],
seq
[
0
],
seq
[
1
]])
seq
=
numpy
.
array
([
seq
[
2
],
seq
[
0
],
seq
[
1
]])
elif
repr
==
'Il'
:
seq
=
N
.
seq
[
2
::
-
1
]
seq
=
seq
[
2
::
-
1
]
elif
repr
==
'IIl'
:
seq
[
1
:
3
]
=
N
.
array
([
seq
[
2
],
seq
[
1
]])
seq
[
1
:
3
]
=
numpy
.
array
([
seq
[
2
],
seq
[
1
]])
elif
repr
==
'IIIl'
:
seq
[
0
:
2
]
=
N
.
array
([
seq
[
1
],
seq
[
0
]])
seq
[
0
:
2
]
=
numpy
.
array
([
seq
[
1
],
seq
[
0
]])
elif
repr
!=
'IIIr'
:
print
'unknown representation'
diag
.
array
=
N
.
take
(
diag
.
array
,
seq
)
diag
.
array
=
numpy
.
take
(
diag
.
array
,
seq
)
return
Transformation
.
Rotation
(
diag
)
*
Transformation
.
Translation
(
-
cm
)
def
applyTransformation
(
self
,
t
):
...
...
@@ -578,7 +579,7 @@ class GroupOfAtoms(object):
universe
=
self
.
universe
()
if
universe
is
None
:
raise
ValueError
(
"object not in a universe"
)
array
=
N
.
zeros
((
universe
.
numberOfAtoms
(),),
N
.
I
nt
)
array
=
numpy
.
zeros
((
universe
.
numberOfAtoms
(),),
numpy
.
i
nt
)
mask
=
ParticleProperties
.
ParticleScalar
(
universe
,
array
)
for
a
in
self
.
atomIterator
():
mask
[
a
]
=
1
...
...
@@ -708,14 +709,14 @@ class Collection(GroupOfAtoms, Visualization.Viewable):
within the rectangular volume
:rtype: :class:`~MMTK.Collections.Collection`
"""
x1
=
N
.
minimum
(
p1
.
array
,
p2
.
array
)
x2
=
N
.
maximum
(
p1
.
array
,
p2
.
array
)
x1
=
numpy
.
minimum
(
p1
.
array
,
p2
.
array
)
x2
=
numpy
.
maximum
(
p1
.
array
,
p2
.
array
)
in_box
=
[]
for
o
in
self
.
objects
:
r
=
o
.
position
().
array
if
N
.
logical_and
.
reduce
(
\
N
.
logical_and
(
N
.
less_equal
(
x1
,
r
),
N
.
less
(
r
,
x2
))):
if
numpy
.
logical_and
.
reduce
(
\
numpy
.
logical_and
(
numpy
.
less_equal
(
x1
,
r
),
numpy
.
less
(
r
,
x2
))):
in_box
.
append
(
o
)
return
Collection
(
in_box
)
...
...
@@ -937,9 +938,9 @@ class PartitionedCollection(Collection):
self
.
all
=
None
def
partitionIndex
(
self
,
x
):
return
(
int
(
N
.
floor
(
x
[
0
]
/
self
.
partition_size
)),
int
(
N
.
floor
(
x
[
1
]
/
self
.
partition_size
)),
int
(
N
.
floor
(
x
[
2
]
/
self
.
partition_size
)))
return
(
int
(
numpy
.
floor
(
x
[
0
]
/
self
.
partition_size
)),
int
(
numpy
.
floor
(
x
[
1
]
/
self
.
partition_size
)),
int
(
numpy
.
floor
(
x
[
2
]
/
self
.
partition_size
)))
def
objectList
(
self
):
return
sum
(
self
.
partition
.
values
(),
[
self
.
undefined
])
...
...
@@ -979,7 +980,7 @@ class PartitionedCollection(Collection):
y
=
int
(
round
(
point
[
1
]
/
self
.
partition_size
))
z
=
int
(
round
(
point
[
2
]
/
self
.
partition_size
))
d
=
(
Vector
(
x
,
y
,
z
)
*
self
.
partition_size
-
point
).
length
()
n
=
int
(
N
.
ceil
((
edge
+
d
)
/
(
2.
*
self
.
partition_size
)))
n
=
int
(
numpy
.
ceil
((
edge
+
d
)
/
(
2.
*
self
.
partition_size
)))
objects
=
[]
for
nx
in
range
(
-
n
,
n
):
for
ny
in
range
(
-
n
,
n
):
...
...
@@ -996,11 +997,11 @@ class PartitionedCollection(Collection):
minsq
=
min
**
2
maxsq
=
max
**
2
for
index
in
self
.
partition
.
keys
():
d1
=
self
.
partition_size
*
N
.
array
(
index
)
-
point
.
array
d1
=
self
.
partition_size
*
numpy
.
array
(
index
)
-
point
.
array
d2
=
d1
+
self
.
partition_size
dmin
=
(
d1
>
0.
)
*
d1
-
(
d2
<
0.
)
*
d2
dminsq
=
N
.
add
.
reduce
(
dmin
**
2
)
dmaxsq
=
N
.
add
.
reduce
(
N
.
maximum
(
d1
**
2
,
d2
**
2
))
dminsq
=
numpy
.
add
.
reduce
(
dmin
**
2
)
dmaxsq
=
numpy
.
add
.
reduce
(
numpy
.
maximum
(
d1
**
2
,
d2
**
2
))
if
dminsq
>=
minsq
and
dmaxsq
<=
maxsq
:
objects
.
addObject
(
self
.
partition
[
index
])
elif
dmaxsq
>=
minsq
and
dminsq
<=
maxsq
:
...
...
@@ -1025,17 +1026,17 @@ class PartitionedCollection(Collection):
for
o1
,
o2
in
Utility
.
pairs
(
zip
(
objects
,
pos
)):
if
(
o2
[
1
]
-
o1
[
1
]).
length
()
<=
cutoff
:
pairs
.
append
((
o1
[
0
],
o2
[
0
]))
partition_cutoff
=
int
(
N
.
floor
((
cutoff
/
self
.
partition_size
)
**
2
))
ones
=
N
.
array
([
1
,
1
,
1
])
zeros
=
N
.
array
([
0
,
0
,
0
])
partition_cutoff
=
int
(
numpy
.
floor
((
cutoff
/
self
.
partition_size
)
**
2
))
ones
=
numpy
.
array
([
1
,
1
,
1
])
zeros
=
numpy
.
array
([
0
,
0
,
0
])
keys
=
self
.
partition
.
keys
()
for
i
in
range
(
len
(
keys
)):
p1
=
keys
[
i
]
for
j
in
range
(
i
+
1
,
len
(
keys
)):
p2
=
keys
[
j
]
d
=
N
.
maximum
(
abs
(
N
.
array
(
p2
)
-
N
.
array
(
p1
))
-
d
=
numpy
.
maximum
(
abs
(
numpy
.
array
(
p2
)
-
numpy
.
array
(
p1
))
-
ones
,
zeros
)
if
N
.
add
.
reduce
(
d
*
d
)
<=
partition_cutoff
:
if
numpy
.
add
.
reduce
(
d
*
d
)
<=
partition_cutoff
:
for
o1
,
pos1
in
zip
(
self
.
partition
[
p1
],
positions
[
p1
]):
for
o2
,
pos2
in
zip
(
self
.
partition
[
p2
],
...
...
MMTK/ConfigIO.py
View file @
3e132477
...
...
@@ -13,12 +13,13 @@ PDB files are handled in :class:`~MMTK.PDB`.
__docformat__
=
'restructuredtext'
import
numpy
from
MMTK
import
PDB
,
Units
,
Utility
from
Scientific.Geometry.Objects3D
import
Sphere
,
Cone
,
Plane
,
Line
,
\
rotatePoint
from
Scientific.Geometry
import
Vector
from
Scientific.Visualization
import
VRML
from
Scientific
import
N
import
os
#
...
...
@@ -164,8 +165,8 @@ class VRMLWireframeFile(VRML.VRMLFile):
self
.
warning
=
0
self
.
color_values
=
color_values
if
self
.
color_values
is
not
None
:
lower
=
N
.
minimum
.
reduce
(
color_values
.
array
)
upper
=
N
.
maximum
.
reduce
(
color_values
.
array
)
lower
=
numpy
.
minimum
.
reduce
(
color_values
.
array
)
upper
=
numpy
.
maximum
.
reduce
(
color_values
.
array
)
self
.
color_scale
=
VRML
.
ColorScale
((
lower
,
upper
))
def
write
(
self
,
object
,
configuration
=
None
,
distance
=
None
):
...
...
MMTK/DCD.py
View file @
3e132477
...
...
@@ -17,10 +17,10 @@ format and is therefore platform-dependent.
__docformat__
=
'restructuredtext'
import
numpy
import
MMTK_DCD
from
MMTK
import
PDB
,
Trajectory
,
Units
from
Scientific
import
N
class
DCDReader
(
Trajectory
.
TrajectoryGenerator
):
...
...
@@ -73,9 +73,9 @@ def writeDCD(vector_list, dcd_file_name, factor, atom_order=None,
universe
=
vector_list
[
0
].
universe
natoms
=
universe
.
numberOfPoints
()
if
atom_order
is
None
:
atom_order
=
N
.
array
range
(
natoms
)
atom_order
=
numpy
.
a
range
(
natoms
)
else
:
atom_order
=
N
.
array
(
atom_order
)
atom_order
=
numpy
.
array
(
atom_order
)
i_start
=
0
# always start at frame 0
n_savc
=
1
# save every frame
...
...
@@ -85,9 +85,9 @@ def writeDCD(vector_list, dcd_file_name, factor, atom_order=None,
if
conf_flag
:
vector
=
universe
.
contiguousObjectConfiguration
(
None
,
vector
)
array
=
factor
*
vector
.
array
x
=
N
.
take
(
array
[:,
0
],
atom_order
).
astype
(
N
.
F
loat16
)
y
=
N
.
take
(
array
[:,
1
],
atom_order
).
astype
(
N
.
F
loat16
)
z
=
N
.
take
(
array
[:,
2
],
atom_order
).
astype
(
N
.
F
loat16
)
x
=
numpy
.
take
(
array
[:,
0
],
atom_order
).
astype
(
numpy
.
f
loat16
)
y
=
numpy
.
take
(
array
[:,
1
],
atom_order
).
astype
(
numpy
.
f
loat16
)
z
=
numpy
.
take
(
array
[:,
2
],
atom_order
).
astype
(
numpy
.
f
loat16
)
MMTK_DCD
.
writeDCDStep
(
fd
,
x
,
y
,
z
)
MMTK_DCD
.
writeCloseDCD
(
fd
)
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment