Source code for wrapper_autofrag
#
# @BEGIN LICENSE
#
# Psi4: an open-source quantum chemistry software package
#
# Copyright (c) 2007-2016 The Psi4 Developers.
#
# The copyrights for code used from other parties are included in
# the corresponding files.
#
# This program 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 2 of the License, or
# (at your option) any later version.
#
# This program 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, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# @END LICENSE
#
from __future__ import print_function
from __future__ import absolute_import
import math
import psi4
def _autofragment_convert(p, symbol):
# Finding radii for auto-fragmenter
if symbol[p] == 'H':
d = 1.001
if symbol[p] == 'He':
d = 1.012
if symbol[p] == 'Li':
d = 0.825
if symbol[p] == 'Be':
d = 1.408
if symbol[p] == 'B':
d = 1.485
if symbol[p] == 'C':
d = 1.452
if symbol[p] == 'N':
d = 1.397
if symbol[p] == 'O':
d = 1.342
if symbol[p] == 'F':
d = 1.287
if symbol[p] == 'Ne':
d = 1.243
if symbol[p] == 'Na':
d = 1.144
if symbol[p] == 'Mg':
d = 1.364
if symbol[p] == 'Al':
d = 1.639
if symbol[p] == 'Si':
d = 1.716
if symbol[p] == 'P':
d = 1.705
if symbol[p] == 'S':
d = 1.683
if symbol[p] == 'Cl':
d = 1.639
if symbol[p] == 'Ar':
d = 1.595
return d / 1.5
[docs]def auto_fragments(**kwargs):
r"""Detects fragments if the user does not supply them.
Currently only used for the WebMO implementation of SAPT.
:returns: :ref:`Molecule<sec:psimod_Molecule>`) |w--w| fragmented molecule.
:type molecule: :ref:`molecule <op_py_molecule>`
:param molecule: ``h2o`` || etc.
The target molecule, if not the last molecule defined.
:examples:
>>> # [1] replicates with cbs() the simple model chemistry scf/cc-pVDZ: set basis cc-pVDZ energy('scf')
>>> molecule mol {\nH 0.0 0.0 0.0\nH 2.0 0.0 0.0\nF 0.0 1.0 0.0\nF 2.0 1.0 0.0\n}
>>> print mol.nfragments() # 1
>>> fragmol = auto_fragments()
>>> print fragmol.nfragments() # 2
"""
# Make sure the molecule the user provided is the active one
molecule = kwargs.pop('molecule', psi4.get_active_molecule())
molecule.update_geometry()
molname = molecule.name()
geom = molecule.save_string_xyz()
numatoms = molecule.natom()
VdW = [1.2, 1.7, 1.5, 1.55, 1.52, 1.9, 1.85, 1.8]
symbol = list(range(numatoms))
X = [0.0] * numatoms
Y = [0.0] * numatoms
Z = [0.0] * numatoms
Queue = []
White = []
Black = []
F = geom.split('\n')
for f in range(numatoms):
A = F[f + 1].split()
symbol[f] = A[0]
X[f] = float(A[1])
Y[f] = float(A[2])
Z[f] = float(A[3])
White.append(f)
Fragment = [[] for i in range(numatoms)] # stores fragments
start = 0 # starts with the first atom in the list
Queue.append(start)
White.remove(start)
frag = 0
while((len(White) > 0) or (len(Queue) > 0)): # Iterates to the next fragment
while(len(Queue) > 0): # BFS within a fragment
for u in Queue: # find all nearest Neighbors
# (still coloured white) to vertex u
for i in White:
Distance = math.sqrt((X[i] - X[u]) * (X[i] - X[u]) +
(Y[i] - Y[u]) * (Y[i] - Y[u]) +
(Z[i] - Z[u]) * (Z[i] - Z[u]))
if Distance < _autofragment_convert(u, symbol) + _autofragment_convert(i, symbol):
Queue.append(i) # if you find you, put it in the que
White.remove(i) # and remove it from the untouched list
Queue.remove(u) # remove focus from Queue
Black.append(u)
Fragment[frag].append(int(u + 1)) # add to group (adding 1 to start
# list at one instead of zero)
if(len(White) != 0): # cant move White->Queue if no more exist
Queue.append(White[0])
White.remove(White[0])
frag += 1
new_geom = """\n"""
for i in Fragment[0]:
new_geom = new_geom + F[i].lstrip() + """\n"""
new_geom = new_geom + """--\n"""
for j in Fragment[1]:
new_geom = new_geom + F[j].lstrip() + """\n"""
new_geom = new_geom + """units angstrom\n"""
moleculenew = psi4.Molecule.create_molecule_from_string(new_geom)
moleculenew.set_name(molname)
moleculenew.update_geometry()
moleculenew.print_cluster()
psi4.print_out(""" Exiting auto_fragments\n""")
return moleculenew