From 319617dafc3055a5f01ceab04ece85c2761e15b2 Mon Sep 17 00:00:00 2001 From: trisyoungs Date: Thu, 5 Dec 2013 14:18:15 +0000 Subject: [PATCH] Modified Pattern::createMatrices() to construct only minimal connectivity matrix for larger patterns (>999 atoms) since the (no doubt sub-optimal) method to calculate the full connectiviy matrix involves a triple loop over atoms. Rewrote part of ato export filter to be far more efficient, and thus not fail for large systems, and added option to omit rotational headgroup output. --- CMakeLists.txt | 2 +- aten.spec | 2 +- configure.ac | 2 +- data/filters/ato | 131 ++++++++++++++++++++++++---------------- extra/aten.dsc | 4 +- src/base/pattern.cpp | 141 ++++++++++++++++++++++++++++++++----------- src/main/version.h | 4 +- 7 files changed, 191 insertions(+), 95 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b272a09be..06c35fd0c 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ project(Aten) set(DESCRIPTION "Aten - Atomic configuration builder and editor") set(AUTHOR "Tristan Youngs") set(VERSION_MAJOR "1") -set(VERSION_MINOR "857") +set(VERSION_MINOR "858") set(VERSION_PATCH "1") set(CMAKE_BUILD_TYPE "Release") diff --git a/aten.spec b/aten.spec index b1c9e9b78..89a48f672 100644 --- a/aten.spec +++ b/aten.spec @@ -4,7 +4,7 @@ # Name, brief description, and version Summary: Aten - Atomic configuration builder and editor Name: %{shortname} -Version: 1.857 +Version: 1.858 Release: 1 License: GPL %define fullname %{name}-%{version} diff --git a/configure.ac b/configure.ac index c05b4af06..492757dba 100644 --- a/configure.ac +++ b/configure.ac @@ -3,7 +3,7 @@ AC_PREREQ(2.60) # Set program name, version, bug-address and source directory -m4_define([ATEN_VERSION],[1.857]) +m4_define([ATEN_VERSION],[1.858]) AC_INIT(aten,ATEN_VERSION,tris@projectaten.net) AC_CONFIG_SRCDIR([src/main.cpp]) diff --git a/data/filters/ato b/data/filters/ato index 5f4914ed9..b6c041b25 100644 --- a/data/filters/ato +++ b/data/filters/ato @@ -2,6 +2,7 @@ # Created: 23/03/2011 # Last modified: 02/09/2011 # ChangeLog: +# 05/12/2013 - Restraint calculation now much more efficient. Added option to control output of headgroup rotations. # 02/09/2011 - Tweaked to use capitalised command/member names in v1.8 # 28/05/2011 - Added writing of random numbers on export. Fixed import to recognise end of LJ specifications. # 23/03/2011 - Initial version. @@ -134,10 +135,10 @@ filter(type="exportmodel",name="EPSR ATO File", nickname="ato", extension="ato", int nmols, n, m, o, mol, nresb, nresa, nrest, nrot, idj, count; double rij, mass, boxsize; Vector com, r; - Atom i, j; + Atom i, j, k, l; Bound b; - Bond bnd; - String rotText, s; + Bond b1, b2, b3; + String rotText, resText, s; # Main dialog creation function void createDefaultDialog(Dialog ui) @@ -158,6 +159,11 @@ filter(type="exportmodel",name="EPSR ATO File", nickname="ato", extension="ato", group.addDoubleSpin("epsr_mtrans", "Whole Molecule Translations", 0.0, 10.0, 0.01, 0.1); ui.addDoubleSpin("epsr_vtemp", "Vibrational Temperature", 0.0, 1500.0, 10.0, 65.0); + # Output Options + group = ui.addGroup("outputoptions", "Output Options", -1, -1, 1); + group.verticalFill = TRUE; + group.addCheck("outputrotations", "Write Rotational Groups", 1); + # Restraints group = ui.addGroup("restraints", "Restraints", -1, -1, 1); ui.addRadioGroup("restraintgroup"); @@ -168,8 +174,8 @@ filter(type="exportmodel",name="EPSR ATO File", nickname="ato", extension="ato", ui.addRadioGroup("connectivitygroup"); w.verticalFill = TRUE; w.addRadioButton("restrain_bonds", "Across Bonds", "connectivitygroup", 1); - w.addRadioButton("restrain_angles", "Across Angles","connectivitygroup", 0); - w.addRadioButton("restrain_torsions", "Across Torsions", "connectivitygroup", 0); + w.addRadioButton("restrain_angles", "Across Bonds and Angles","connectivitygroup", 0); + w.addRadioButton("restrain_torsions", "Across Bonds, Angles and Torsions", "connectivitygroup", 0); w.addDoubleSpin("restrain_con_max", "Maximum Distance", 0.0, 10.0, 0.1, 2.4); w = group.addFrame("restrain_dist_frame", 2, 2); @@ -184,6 +190,12 @@ filter(type="exportmodel",name="EPSR ATO File", nickname="ato", extension="ato", if (!showDefaultDialog()) error("Canceled through dialog.\n"); Dialog ui = defaultDialog(); + // Grab some dialog values + int resType = ui.asInteger("restrain_bonds"); + if (ui.asInteger("restrain_angles")) resType = 2; + else if (ui.asInteger("restrain_torsions")) resType = 3; + double rijmax = ui.asDouble("restrain_con_max"); + // First, some checks. We need valid patterns and a cubic cell. if (!createExpression(FALSE,TRUE)) error("Error: Can't write ATO file without valid forcefield types assigned to all atoms.\n"); @@ -250,58 +262,71 @@ filter(type="exportmodel",name="EPSR ATO File", nickname="ato", extension="ato", // Restraint information count = 0; + resText = ""; if (ui.asInteger("restrain_con")) { - nresb = 0; - nresa = 0; - nrest = 0; - double rijmax = ui.asDouble("restrain_con_max"); - if (ui.asInteger("restrain_bonds")) for (b in p.bonds) if (((b.id[1] == n) || (b.id[2] == n)) && (geometry(p.atoms[b.id[1]],p.atoms[b.id[2]]) < rijmax)) ++nresb; - if (ui.asInteger("restrain_angles")) for (b in p.angles) if (((b.id[1] == n) || (b.id[3] == n)) && (geometry(p.atoms[b.id[1]],p.atoms[b.id[3]]) < rijmax)) ++nresa; - if (ui.asInteger("restrain_torsions")) for (b in p.torsions) if (((b.id[1] == n) || (b.id[4] == n)) && (geometry(p.atoms[b.id[1]],p.atoms[b.id[4]]) < rijmax)) ++nrest; - writeLineF(" %2i ", nresa+nresb+nrest); - if (nresb > 0) for (b in p.bonds) + // Loop over bonds to this atom + for (b1 in i.bonds) { - // Search for bonds in which this atom is involved - if (b.id[1] == n) idj = b.id[2]; - else if (b.id[2] == n) idj = b.id[1]; - else continue; - rij = geometry(p.atoms[b.id[1]],p.atoms[b.id[2]]); - if (rij > rijmax) continue; - if (count && (count%5 == 0)) writeLineF("\n"); - ++count; - writeLineF("%4i %9.3e ", idj, rij); + j = b1.partner(i); + + // Always restrain bond distance (if within max distance allowed).... + rij = geometry(b1.i.id, b1.j.id); + if (rij < rijmax) + { + ++count; + sprintf(s, "%4i %9.3e ", j.id - p.firstAtomId + 1, rij); + resText += s; + if (count%5 == 0) resText += "\n"; + } + + // Restrain across angles? + if (resType > 1) + { + // Loop over bonds on j + for (b2 in j.bonds) + { + if (b1 == b2) continue; + k = b2.partner(j); + rij = geometry(k.id, i.id); + if (rij < rijmax) + { + ++count; + sprintf(s, "%4i %9.3e ", k.id - p.firstAtomId + 1, rij); + resText += s; + if (count%5 == 0) resText += "\n"; + } + + // Restrain across torsions? + if (resType > 2) + { + // Loop over bonds on k + for (b3 in k.bonds) + { + if (b2 == b3) continue; + l = b3.partner(k); + rij = geometry(l.id, i.id); + if (rij < rijmax) + { + ++count; + sprintf(s, "%4i %9.3e ", l.id - p.firstAtomId + 1, rij); + resText += s; + if (count%5 == 0) resText += "\n"; + } + + } + } + } + } + + } - if (nresa > 0) for (b in p.angles) - { - // Search for angles in which this atom is involved - if (b.id[1] == n) idj = b.id[3]; - else if (b.id[3] == n) idj = b.id[1]; - else continue; - rij = geometry(p.atoms[b.id[1]],p.atoms[b.id[3]]); - if (rij > rijmax) continue; - if (count && (count%5 == 0)) writeLineF("\n"); - ++count; - writeLineF("%4i %9.3e ", idj, rij); - } - if (nrest > 0) for (b in p.torsions) - { - // Search for torsions in which this atom is involved - if (b.id[1] == n) idj = b.id[4]; - else if (b.id[4] == n) idj = b.id[1]; - else continue; - rij = geometry(p.atoms[b.id[1]],p.atoms[b.id[4]]); - if (rij > rijmax) continue; - if (count && (count%5 == 0)) writeLineF("\n"); - ++count; - writeLineF("%4i %9.3e ", idj, rij); - } - writeLineF("\n"); + + writeLineF(" %2i %s\n", count, resText); } else { nresb = 0; - double rijmax = ui.asDouble("restrain_dist_max"); for (idj = 1; idj<=p.nMolAtoms; ++idj) if ((n != idj) && (geometry(p.atoms[n],p.atoms[idj]) < rijmax)) ++nresb; writeLineF(" %4i", nresb); if (nresb > 0) for (idj = 1; idj<=p.nMolAtoms; ++idj) @@ -320,14 +345,14 @@ filter(type="exportmodel",name="EPSR ATO File", nickname="ato", extension="ato", // We will construct a list of rotations and count them up as we go. A string will be created containing all the info to be written out. nrot = 0; rotText = ""; - for (b in p.bonds) + if (ui.asInteger("outputrotations")) for (b in p.bonds) { if (p.atomsInRing(b.id[1],b.id[2])) continue; if (p.atoms[b.id[1]].nBonds == 1) continue; if (p.atoms[b.id[2]].nBonds == 1) continue; - bnd = p.atoms[b.id[1]].findBond(p.atoms[b.id[2]]); + b1 = p.atoms[b.id[1]].findBond(p.atoms[b.id[2]]); srcmodel.selectNone(); - srcmodel.selectTree(p.atoms[b.id[1]], bnd); + srcmodel.selectTree(p.atoms[b.id[1]], b1); // Check here the number of selected atoms - if greater than half the atoms in the molecule then we're better off with the inverse selection! if (srcmodel.nSelected > 0.5*p.nMolAtoms) { @@ -399,7 +424,7 @@ filter(type="exportmodel",name="EPSR ATO File", nickname="ato", extension="ato", // Extra data // Used by fmole to keep non-bonded atoms apart - writeLineF(" %10.4e %10.4e\n", 1.0, 3.0); + writeLineF(" %10.4e %10.4e\n", 1.0, 1.0); // Random numbers for restart purposes for (n=0; n<15; ++n) writeLineF(" %i", randomI()); diff --git a/extra/aten.dsc b/extra/aten.dsc index 6d505441b..2601c9175 100644 --- a/extra/aten.dsc +++ b/extra/aten.dsc @@ -1,9 +1,9 @@ Format: 1.0 Source: aten -Version: 1.857 +Version: 1.858 Binary: aten Maintainer: Tristan Youngs Architecture: any Build-Depends: debhelper (>= 4.1.16), libqt4-dev | libqt4-core, libqt4-opengl-dev, libreadline5-dev | libreadline-dev, libgl1-mesa-dev, pkgconfig | pkg-config, libncurses5 Files: - 0f7cda0fc1d5e6e2ca1c6247eb944c61 4309758 aten-1.857.tar.gz + 0f7cda0fc1d5e6e2ca1c6247eb944c61 4309758 aten-1.858.tar.gz diff --git a/src/base/pattern.cpp b/src/base/pattern.cpp index 494080da2..22885071a 100644 --- a/src/base/pattern.cpp +++ b/src/base/pattern.cpp @@ -722,16 +722,8 @@ void Pattern::createMatrices() for (n=0; nnext) - { - conMatrix_[ pb->atomId(0) ] [ pb->atomId(1) ] = 1; - conMatrix_[ pb->atomId(1) ] [ pb->atomId(0) ] = 1; - } + for (n=0; nnext) { - if (conMatrix_[a1][a2] != 0) + conMatrix_[ pb->atomId(0) ] [ pb->atomId(1) ] = 1; + conMatrix_[ pb->atomId(1) ] [ pb->atomId(0) ] = 1; + } + + msg.print("transforming (full)....."); + // Now, transform into the connectivity matrix. + for (a1=0; a1a1 and a1->a2 is less than m->a2 - // The last check means that only the minimum distance m->a2 persists at the end - if ((conMatrix_[m][a2] == 0) || (conMatrix_[a1][m]+conMatrix_[a1][a2] < conMatrix_[m][a2])) - conMatrix_[m][a2] = conMatrix_[a1][m] + conMatrix_[a1][a2]; + // We only potentially increase the distance if : + // 1) The atom 'm' is *not equal* to a2 (i.e. is not itself) + // 2) The atom 'm' we're looking at is also bound to a1. + if ((m != a2) && (conMatrix_[a1][m] != 0)) + if ((conMatrix_[a1][m] != 0)) + { + // We only actually increase the distance if : + // 1) Atom 'm' is not directly bound to a2 **OR** + // 2) The combined distances of m->a1 and a1->a2 is less than m->a2 + // The last check means that only the minimum distance m->a2 persists at the end + if ((conMatrix_[m][a2] == 0) || (conMatrix_[a1][m]+conMatrix_[a1][a2] < conMatrix_[m][a2])) + conMatrix_[m][a2] = conMatrix_[a1][m] + conMatrix_[a1][a2]; + } } } } } } + else + { + // Create minimal transformation matrix, using only bond, angle, and torsion data + msg.print("transforming (minimal)....."); + + // There may be more than one consecutive bound fragment in the pattern, so we must perform treeSelects in order to populate the initial matrix + Atom* i = firstAtom_; + int count = 100, ii, jj, diagii; + while (i != NULL) + { + // Treeselect from current atom + parent_->selectTree(i, TRUE); + + // For the current marked selection, set the diagonal matrix elements to the current 'count' value + for (Refitem* ri = parent_->selection(TRUE); ri != NULL; ri = ri->next) + { + ii = ri->item->id() - startAtom_; + conMatrix_[ii][ii] = count; + } + + // Find next unmarked atom + while (i && i->isSelected(TRUE)) i = i->next; + + // Deselect all atoms, and increase count + parent_->selectNone(TRUE); + ++count; + } + + // The diagonal elements of conMatrix_ now indicate the parent fragments of each atom + // For all atoms within a given fragment, set the connectivity to the diagonal value to start with + for (ii = 0; ii < nAtoms_; ++ii) + { + diagii = conMatrix_[ii][ii]; + for (jj = 0; jj < nAtoms_; ++jj) + { + if (diagii != conMatrix_[jj][jj]) continue; + conMatrix_[ii][jj] = diagii; + conMatrix_[jj][ii] = diagii; + } + } + + // Done with the diagonals now, so zero them + for (ii = 0; ii < nAtoms_; ++ii) conMatrix_[ii][ii] = 0; + + // Now, add bonds to matrix + for (pb = bonds_.first(); pb != NULL; pb = pb->next) + { + conMatrix_[ pb->atomId(0) ] [ pb->atomId(1) ] = 1; + conMatrix_[ pb->atomId(1) ] [ pb->atomId(0) ] = 1; + } + + // Angles (but don't overwrite bonds) + for (pb = angles_.first(); pb != NULL; pb = pb->next) + { + if (conMatrix_[pb->atomId(0)][pb->atomId(2)] == 1) continue; + conMatrix_[pb->atomId(0)][pb->atomId(2)] = 2; + conMatrix_[pb->atomId(2)][pb->atomId(0)] = 2; + } + + // Torsions (but don't overwrite angles or bonds) + for (pb = torsions_.first(); pb != NULL; pb = pb->next) + { + if (conMatrix_[pb->atomId(0)][pb->atomId(3)] < 3) continue; + conMatrix_[pb->atomId(0)][pb->atomId(3)] = 3; + conMatrix_[pb->atomId(3)][pb->atomId(0)] = 3; + } + } msg.print("done.\n"); -/* printf("Connectivity Matrix\n"); - for (n=0; n