Introduction
This is a quick guide in evolution of "KICK" stochastic search procedure.
It has had many iterations in the past 16 years in our research group (Vrček-Šakić), including moving from FORTRAN to bash for usage in command line without compiler, and again rewritten in javascript for usage in web-pages.
Two other methods briefly described here are:
Coalescence Kick in c++ by Boris Averkiev and
SnippetKick in Perl by Osvaldo Yañez Osses.
The original paper was published in 2004, and currently (2020) it has been cited 177 times.
"Stochastic Search for Isomers on a Quantum Mechanical Surface"
M. Saunders, J. Comput. Chem. 2004, 25, 621 - 626.
Main procedure is perfectly summed in prof. Saunders own words:
The procedure can begin with any structure for the molecule. It is submitted to GAUSSIAN optimization. The minimum energy structure obtained is then subject to an operation called a kick. Each atom is moved a random distance in a random direction. The only limitation is the maximum distance the atoms are to be moved. Each atom is kicked to a position within a sphere of radius R around its initial position, where R is the maximum kick distance. Calculating such randomly moved positions is extremely simple. One generates random values between -R and +R for the displacements along the x, y, and z axes. An atom kicked in such a way will, in fact, wind up within a cube 2R on a side instead of a sphere of radius R. The orientation of this cube would be changed by rotating the molecule in Cartesian space. To have kicked positions fall within a sphere instead of a cube, one need only calculate the distance moved [sqrt(x^2 + y^2 - z^2)]. If it is greater than R, this kick is rejected and a new set of three random displacements is chosen. Effectively, we are cutting off the corners of the cube and making it into a sphere.
Here an attempt to recreate minimal steps needed for javascript implementation of the method will be undertaken with practical code snippets that can be useful in other projects. Additionally, different variations routed in the original procedure will be discussed.
Original Saunders KICK version
Directly in javascript this would look something like this:
//definition of parameters
var param = 2.0 //maximum distance parameter R
var checksphere = 0 //sphere test
//start of while checksphere loop
while (checkspehe == 0) {
var dx = (Math.random()*2 - 1)*Math.random();
var dy = (Math.random()*2 - 1)*Math.random();
var dz = (Math.random()*2 - 1)*Math.random();
var size = Math.sqrt(Math.pow(dx,2)+Math.pow(dy,2)+Math.pow(dz,2));
//checking if vector size is less then maximum distance parameter R
if (size <= param) {
checksphere = 1; //exit from the loop if true
} else {
checksphere = 0; //redo vector generation
}
//finishing checksphere loop
}
Although this is a good way to construct such a thing, this approach is more mathematically sound:
It consist of construction of a unit (direction) vector which is generated from randomly chosen
values for x, y, and z coordinates in the range of -1 to +1, and then normalized, by calculating
vector length, and then dividing each coordinate by said amount. Loop is avoided if things are
done this way. In javascript, current procedure looks as follows:
var param = 2.0 //maximum distance parameter R
var minDist = 0.1 //minimal distance parameter
//to find a random x part of the vector (vx), first part a random +/- sign via
//(Math.round(Math.random()) * 2 - 1) which multiplies random number generated by
//Math.random(), and additional 0.0001 just to avoid getting rounding errors.
var vx = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vy = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vz = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
//generated radnom numbers are squared
var vx2 = Math.pow(vx,2);
var vy2 = Math.pow(vy,2);
var vz2 = Math.pow(vz,2);
//vector length from vx, vy, and vz is calculated
var vsum = vx2 + vy2 + vz2;
var vecnorm = Math.sqrt(vsum);
var vecsize = Math.random()*(param) + minDist;
var vectorX = (vx / vecnorm) * vecsize ;
var vectorY = (vy / vecnorm) * vecsize ;
var vectorZ = (vz / vecnorm) * vecsize ;
"After all the atoms are randomly moved in this way, quantum mechanical optimization is carried out again." Movement of atoms here is straightforward. After coordinates are read from text-input, or from a file, a loop is made to generate random vector for each atom, and new coordinates are saved in an array. If all checks are successful, that array is written to output, either as text in a web-page and/or in a file for download. If our input area was converted to an array called texts, we can extract coordinates in .xyz style in this way:
//definition of variables
var coorA = [];
var coorArrX = [];
var coorArrY = [];
var coorArrZ = [];
//reading from texts variable
for (var i=0; i < texts.length; i++) {
var coor = texts[i].replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").split(" ");
//populating variables
coorA[i] = coor[0];
coorArrX[i] = Number(coor[1]);
coorArrY[i] = Number(coor[2]);
coorArrZ[i] = Number(coor[3]);
}
//parameters needed for calculation
var param = 2.0 //maximum distance parameter R
var minDist = 0.1 //minimal distance parameter
//definition of new variables
var newA = [];
var newX = [];
var newY = [];
var newZ = [];
Random generation of vector and addition to it can now be achieved either like this:
//starting atom loop
for (var i=0; i < texts.length; i++) {
var checksphere = 0 //sphere test
while (checkspehe == 0) {
var dx = (Math.random()*2 - 1)*Math.random();
var dy = (Math.random()*2 - 1)*Math.random();
var dz = (Math.random()*2 - 1)*Math.random();
var size = Math.sqrt(Math.pow(dx,2)+Math.pow(dy,2)+Math.pow(dz,2));
//checking if vector size is less then maximum distance parameter R
if (size <= param) {
checksphere = 1; //exit from the loop if true
newA[i] = coorA[i];
newX[i] = coorArrX[i]+dx;
newY[i] = coorArrY[i]+dy;
newZ[i] = coorArrZ[i]+dz;
} else {
checksphere = 0; //redo vector generation
}
//finishing checksphere loop
}
//finishing atom loop
}
Or like this, if different vector generation procedure is used.
//starting atom loop
for (var i=0; i < texts.length; i++) {
var vx = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vy = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vz = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vx2 = Math.pow(vx,2);
var vy2 = Math.pow(vy,2);
var vz2 = Math.pow(vz,2);
var vsum = vx2 + vy2 + vz2;
var vecnorm = Math.sqrt(vsum);
var vecsize = Math.random()*(param) + minDist;
var vectorX = (vx / vecnorm) * vecsize ;
var vectorY = (vy / vecnorm) * vecsize ;
var vectorZ = (vz / vecnorm) * vecsize ;
var dx = vectorX;
var dy = vectorY;
var dz = vectorZ;
newA[i] = coorA[i];
newX[i] = coorArrX[i]+dx;
newY[i] = coorArrY[i]+dy;
newZ[i] = coorArrZ[i]+dz;
//finishing atom loop
}
The first real modification that one has to do to get meaningful starting structures is to do a quick check of inter-atom distances. Here
minDist
var status = []; //status of calculation
var dist = []; //distance variable
//first atom loop
for (var i=0; i < texts.length; i++) {
//second atom loop
for (var k=i+1; k < texts.length; k++) {
//distance calculation (same formula for vector size calculation)
var x2 = Math.pow((newX[i]-newX[k]),2);
var y2 = Math.pow((newY[i]-newY[k]),2);
var z2 = Math.pow((newZ[i]-newZ[k]),2);
var sum = x2 + y2 + z2;
var somedist = Math.sqrt(sum);
//puting somedist in dist array
dist.push(somedist);
//end first atom loop
}
//end second atom loop
}
//finding minimal distance in dist array
var mininal_dist_found = Math.min.apply(Math, dist);
//decistion if criteria is fullfiled
if (mininal_dist_found < minDist) {
status = "BAD";
} else {
status = "OK";
}
Of course, this is only for one instance/try of the complete procedure, and fromsum of covalent radius this procedure a number of starting points for further calculation are expected.
To do this, a larger loop is constructed:
//Here goes input (.xyz coordinates) reading
counter = 10; //generate 10 structures
errorCounter = 10; //error counter
//parameters needed for calculation
var param = 2.0 //maximum distance parameter R
var minDist = 0.1 //minimal distance parameter
//counter loop
for (var count=0; count < counter; count++) {
//adding error variable and how it should change parameter R
if (error1 == errorCounter ) {
error1 = 0;
param = param + 0.1*param;
//here goes vector generation and distance calculation
//errors and decision if this is OK structure goes into this:
//decision if criteria has been fulfilled
if (mininal_dist_found < minDist) {
status = "BAD";
error1++: //add a error counter
count=count-1; //rerun the complete step
} else {
status = "OK";
}
//end of counter loop
}
Coalescence KICK version
Coalescence kick procedure was developed by Boris Averkiev in the Boldyrev group at Utah University."Geometry and Electronic Structure of Doped Clusters via the Coalescence Kick Method" B. Averkiev, 2009, Ph.D. dissertation, Utah State University, Logan, Utah, 18-25
Two major drawback of Saunders original procedure as mentioned in the thesis:
1. "Initially generated structure containing too short atom-atom distances."
2. "Fragmented initially generated structure. The fragmentation of structure results in very slow conversion of SCF. Only in rare cases initial fragmented structure is optimized into a non-fragmented one." Those problems become "main obstacle when we tried to calculate medium size clusters with more than 12 atoms." "To avoid above-mentioned drawbacks of the original version of Kick Method, we designed the modified “Coalescence” Kick Method. In this method initially generated random structure is first checked for connectivity, and then if the structure is fragmented, the coalescence procedure is applied to it. Two atoms are considered to be connected if the inter-atomic distance satisfies the sum of their covalent (or van der Waals) radii." Coalescence is in "push(ing) all fragments of structure to the center of mass simultaneously." But to do this, after random generation of coordinates, a distance matrix is made. Then close atoms are group together as fragment, and complete fragment is then moved by certain amount to th recalculated center of the mass of the complete system. As atoms are moved closer, fragments are getting larger, aka coalescing, until all atoms are in the same fragment. This can be implemented into the original Saunders stochastic search procedure. Generated random coordinates of atoms are:
newA[]newX[]newY[]newZ[]newF[]information of each atom, and find the coordinates of center of mass for the complete molecule.
First covalent (or van der Waals) radius are put in
newR[]
var newRad = [];
var newF = [];
var covRadius = {He:0.28, H:0.31, F:0.57, Ne:0.58, O:0.66, N:0.71, C:0.76, B:0.85, Be:0.96, Cl:1.02, S:1.05, Ar:1.06, P:1.07, S:1.11, Kr:1.16, As:1.19, Ge:1.2, Se:1.2, Br:1.2, Al:1.21, Zn:1.22, Ga:1.22, Ni:1.24, Co:1.26, Li:1.28, Cu:1.32, Hg:1.32, Pt:1.36, Au:1.36, Te:1.38, Cr:1.39, Mg:1.39, Pd:1.39, Sn:1.39, Sb:1.39, I:1.39, Xe:1.4, Po:1.4, Mg:1.41, Ir:1.41, Rh:1.42, In:1.42, Cd:1.44, Os:1.44, Ag:1.45, Tl:1.45, Ru:1.46, Pb:1.46, Tc:1.47, Bi:1.48, At:1.5, Ra:1.5, Re:1.51, V:1.53, Mo:1.54, Ti:1.6, W:1.62, Nb:1.64, Na:1.66, Cm:1.69, Sc:1.7, Ta:1.7, Zr:1.75, Hf:1.28, Ca: 1.76, Am:1.8, Yb:1.87, Lu:1.87, Pu:1.87, Er:1.89, Y:1.9, Tm:1.9, Np:1.9, Dy:1.92, Ho:1.92, Tb:1.94, Sr:1.95, Gd:1.96, U:1.96, Sm:1.98, Eu:1.98, Pm:1.99, Pa:2, Nd:2.01, K:2.03, Pr:2.03, Ce:2.04, Th:2.06, La:2.07, Ba:2.15, Ac:2.15, Rb:2.2, Ra:2.21, Cs:2.44, Fr:2.6};
for (i=0; i < texts.length; i++) {
newRad[i]=covRadius[coorA[i]]; //every atom has a covalent radius
newF[i]=i; //every atom is a fragment
}
Instead of previous distance calculation, where only inter-atomic distances are checked versus a constant value
minDistdistances are checked if they fall between 70% of sum of covalent radius (minimal distance) and 110% of sum of covalent radius in order to say a bond was made between two atoms.
If inter-atomic distance falls into this category, it is added to the fragment of the fist atom.
var dist = [];
//first atom loop
for (var i=0; i < texts.length; i++) {
//second atom loop
for (var k=i+1; k < texts.length; k++) {
var x2 = Math.pow((newX[i]-newX[k]),2);
var y2 = Math.pow((newY[i]-newY[k]),2);
var z2 = Math.pow((newZ[i]-newZ[k]),2);
var sum = x2 + y2 + z2;
var somedist = Math.sqrt(sum);
var bondlengthOK = 1.1*(newRad[i]+newRad[k]); //maximum distance for a bond
var bondlengthBAD = 0.7*(newRad[i]+newRad[k]); //minimum distance without overlap
dist.push(somedist);
//check
if (somedist <= bondlengthBAD){
demo += (" SMALL DISTANCE ERROR ");
err++ //error counter
break //exit from the loop
} else {
if (somedist <= bondlengthOK){
newF[k]=newF[i]; //adding atom k to the fragment of the i
}
//check end
}
//second atom loop end
}
//first atom loop end
}
Finding total center of the system, alongside finding centers of all the fragments is needed for construction of displacement vector.
function average(p,c,i,a){return p + (c/a.length)}; //averaging function
var totalcentarx = newX.reduce(average,0); //center of whole system
var totalcentary = newY.reduce(average,0);
var totalcentarz = newZ.reduce(average,0);
var fragX = []; //sum of all coordinates inside fragment
var fragY = [];
var fragZ = [];
for (var m=0; m < texts.length; m++) {
fragX[m] = [];
fragY[m] = [];
fragZ[m] = [];
for (var i=0; i < texts.length; i++) {
if (newF[i] == l){
fragX[m].push(newX[i]);
fragY[m].push(newY[i]);
fragZ[m].push(newZ[i]);
}
}
}
var fragXcentar = [];
var fragYcentar = [];
var fragZcentar = [];
for (var m=0; m < texts.length; m++) {
fragXcentar.push(fragX[m].reduce(average,0)); //center of each fragment
fragYcentar.push(fragY[m].reduce(average,0));
fragZcentar.push(fragZ[m].reduce(average,0));
}
Finally a vector is constructed from center of each fragment towards center of the whole system.
Vector is scaled to a value that is sum of 0.1 a.u. and ration of distance parameter- R and loop steps,
where loop steps are number of repetitions given for the system to coalescence. Other conditions for coalescence is that system consists of only one fragment, meaning all atoms in the system are connected.
var deltaX = [];
var deltaY = [];
var deltaZ = [];
for (var m=0; m < texts.length; m++) {
var vx=totalcentarx-fragXcentar[m]; //vector from fragment center to system center
var vy=totalcentary-fragYcentar[m];
var vz=totalcentarz-fragZcentar[m];
var vx2 = Math.pow(vx,2);
var vy2 = Math.pow(vy,2);
var vz2 = Math.pow(vz,2);
var vsum = vx2 + vy2 + vz2;
var vecnorm = Math.sqrt(vsum);
deltaX[m]=(0.1+param/loopMax)*(vx/vecnorm); //step of 0.1 a.u. + param/loop in total
deltaY[m]=(0.1+param/loopMax)*(vy/vecnorm);
deltaZ[m]=(0.1+param/loopMax)*(vz/vecnorm);
}
for (i=0; i < texts.length; i++) {
for (m=0; m < texts.length; m++) {
if (newF[i] == m) {
newX[i] = newX[i] + deltaX[m]; //addition of vector to coordinate
newY[i] = newY[i] + deltaY[m];
newZ[i] = newZ[i] + deltaZ[m];
}
}
}
Easier way to (re)count the (remaining) fragments:
function countUnique(iterable) {
return new Set(iterable).size;
}
var fragCount = countUnique(newF);
This code can be fitted to all the rest discussed procedures. If covalent radius is exchanged to van der Waals radius
weak inter-atomic(molecular) bonds can be encouraged. This is tempting, but sacrifice to randomness has to be expected.
Frozen fragment KICK version
First modification of Saunders original kick stochastic search procedure was by defining a frozen fragment. Atoms in the frozen fragment simply do not move. Defining frozen fragment is straightforward. Atoms in .xyz chemical format are defined as:[element] [x] [y] [z][F]
var coorA = [];
var coorArrX = [];
var coorArrY = [];
var coorArrZ = [];
var coorArrF = [];
for (var i=0; i < texts.length; i++) {
var coor = texts[i].replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").split(" ");
coor.push(0);
coorA[i] = coor[0];
coorArrX[i] = Number(coor[1]);
coorArrY[i] = Number(coor[2]);
coorArrZ[i] = Number(coor[3]);
coorArrF[i] = (coor[4]);
}
To incorporate this into the original script, "atom loop" needs to be modified with one
ifstatetment:
//starting atom loop
for (var i=0; i < texts.length; i++) {
//frozen fragment is denoted with F in fifth column
if (coorArrF[i] == "F") {
//frozen atoms
newA[i] = coorA[i];
newX[i] = coorArrX[i];
newY[i] = coorArrY[i];
newZ[i] = coorArrZ[i];
} else {
//non-frozen atoms
var vx = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vy = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vz = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vx2 = Math.pow(vx,2);
var vy2 = Math.pow(vy,2);
var vz2 = Math.pow(vz,2);
var vsum = vx2 + vy2 + vz2;
var vecnorm = Math.sqrt(vsum);
var vecsize = Math.random()*(param) + minDist;
var vectorX = (vx / vecnorm) * vecsize ;
var vectorY = (vy / vecnorm) * vecsize ;
var vectorZ = (vz / vecnorm) * vecsize ;
var dx = vectorX;
var dy = vectorY;
var dz = vectorZ;
newA[i] = coorA[i];
newX[i] = coorArrX[i]+dx;
newY[i] = coorArrY[i]+dy;
newZ[i] = coorArrZ[i]+dz;
}
//finishing atom loop
}
WIGGLE version
The next major modification is not to move individual atoms, but to move, by "kick" completemolecular fragments or molecules.
This is not a new idea; several authors have done simmilar modification, most notably in a
program called SnippetKick based on the work of Addicoat et al. and written in Perl by Osvaldo Yañez Osses in 2018.
Quick description of the program from github repository:
"The SnippetKick algorithm combines local with global search for molecules and clusters. An
extension of the Kick program developed by Addicoat et al. is described in which chemically
sensible molecular fragments are used in an automated stochastic search algorithm.(ref 1-2)
Ref1: J. Comput. Chem. 2009, 30 (1), 57–64.
Ref2: J. Phys. Chem. A 2006, 110 (13), 4287-4290.
While before mentioned program works beutifully, it still needs to be installed.
Here described procedure works from a webpage.
With molecular fragments, several challenges are present. First molecular fragments need
to be identified and correctly assigned. Atoms in the frozen fragment are noted as "0", and
non-frozen fragments are marked with number defining them.
var coorA = [];
var coorArrX = [];
var coorArrY = [];
var coorArrZ = [];
var coorArrF = [];
for (var i=0; i < texts.length; i++) {
var coor = texts[i].replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").split(" ");
//major point is to add a fifth column if it was forgotten; that atom is then frozen
coor.push(0);
coorA[i] = coor[0];
coorArrX[i] = Number(coor[1]);
coorArrY[i] = Number(coor[2]);
coorArrZ[i] = Number(coor[3]);
coorArrF[i] = Number(coor[4]);
}
var fragCount = Math.max.apply(Math, coorArrF); //number of fragments
For each (non-frozen) fragment, a random vector is needed.
for (var m=1; m <= fragCount; m++) {
var vx = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vy = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vz = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vx2 = Math.pow(vx,2);
var vy2 = Math.pow(vy,2);
var vz2 = Math.pow(vz,2);
var vsum = vx2 + vy2 + vz2;
var vecnorm = Math.sqrt(vsum);
var vecsize = Math.random()*(param) + minDist/2;
var vectorX = (vx / vecnorm) * vecsize ;
var vectorY = (vy / vecnorm) * vecsize ;
var vectorZ = (vz / vecnorm) * vecsize ;
//VECTOR IN ARRAY OF VECTORS
var dx = vectorX;
var dy = vectorY;
var dz = vectorZ;
fragDx[m] = dx;
fragDy[m] = dy;
fragDz[m] = dz;
}
And atoms belonging to each fragment are than translated by said fragment vector.
for (var i=0; i < texts.length; i++) {
for (var m=0; m <= fragCount; m++) {
if (Number(coorArrF[i]) == m) {
var ada = coorA[i];
var xdx = Number(coorArrX[i]) + fragDx[m];
var ydy = Number(coorArrY[i]) + fragDy[m];
var zdz = Number(coorArrZ[i]) + fragDz[m];
newA.push(ada);
newX.push(xdx);
newY.push(ydy);
newZ.push(zdz);
}
}
}
Fragments are not spherical and symmetrical, so it is not only movement in random direction
that is needed, as was the case with "kicking" atoms. For rotation, first a center of each
fragment needs to be set.
var sumcentarx = [];
var sumcentary = [];
var sumcentarz = [];
var atomsInFrag = [];
for (var i=0; i < texts.length; i++) {
for (var m=0; m <= fragCount; m++) {
sumcentarx.push(0);
sumcentary.push(0);
sumcentarz.push(0);
atomsInFrag.push(0);
if (coorArrF[i] == m) {
sumcentarx[m] = sumcentarx[m]+newX[i];
sumcentary[m] = sumcentary[m]+newY[i];
sumcentarz[m] = sumcentarz[m]+newZ[i];
atomsInFrag[m] = atomsInFrag[m]+1;
}
}
}
var centarx = [];
var centary = [];
var centarz = [];
for (var m=0; m <= fragCount; m++) {
centarx[m] = sumcentarx[m] / atomsInFrag[m] ;
centary[m] = sumcentary[m] / atomsInFrag[m] ;
centarz[m] = sumcentarz[m] / atomsInFrag[m] ;
}
Now that center is found, a complex rotation function can be implemented. It is based on quaternium equations. A Quaternion describes a rotation in 3D space.
The Quaternion is mathematically defined as Q = x*i + y*j + z*k + w, where (i,j,k) are imaginary basis vectors.
(x,y,z) can be seen as a vector related to the axis of rotation, while the real multiplier, w, is related to the amount of rotation.
More information can be found at "Understanding Quaternions" R. Goldman, Graphical Models 2011, 73(2), 21 - 49. DOI: 10.1016/j.gmod.2010.10.004.
The complete uderlying mechanism is devided in several functions, that are then called:
function Vector(optional) {
if (optional) {
this.x = optional.x
this.y = optional.y
this.z = optional.z
this.w = optional.w
} else {
this.x = 0
this.y = 0
this.z = 0
this.w = 0
}
}
Vector.prototype.QuaternionMultiply = function(vectorB) {
var out = new Vector();
out.w = this.w*vectorB.w - this.x*vectorB.x - this.y*vectorB.y - this.z*vectorB.z;
out.x = this.w*vectorB.x + this.x*vectorB.w + this.y*vectorB.z - this.z*vectorB.y;
out.y = this.w*vectorB.y - this.x*vectorB.z + this.y*vectorB.w + this.z*vectorB.x;
out.z = this.w*vectorB.z + this.x*vectorB.y - this.y*vectorB.x + this.z*vectorB.w;
this.x = out.x;
this.y = out.y;
this.z = out.z;
this.w = out.w;
}
Vector.prototype.rotate = function (inputaxis, inputangle)
{
var vector = new Vector(this);
vector.w = 0;
var axis = new Vector({
x: inputaxis.x * Math.sin(inputangle/2),
y: inputaxis.y * Math.sin(inputangle/2),
z: inputaxis.z * Math.sin(inputangle/2),
w: Math.cos(inputangle/2)}
);
var axisInv = new Vector({ x: -axis.x, y: -axis.y, z: -axis.z, w: axis.w} );
axis.QuaternionMultiply(vector);
axis.QuaternionMultiply(axisInv);
this.x = axis.x.toFixed(5);
this.y = axis.y.toFixed(5);
this.z = axis.z.toFixed(5);
}
Vector.prototype.length = function () {
var length = Math.sqrt( this.x*this.x + this.y*this.y + this.z*this.z );
return length;
}
Vector.prototype.scale = function (scale) {
this.x *= scale
this.y *= scale
this.z *= scale
}
Vector.prototype.normalize = function() {
var lengthval = this.length()
if (lengthval != 0) {
this.x /= lengthval;
this.y /= lengthval;
this.z /= lengthval;
return true
} else {
return false
}
}
When all functions are declared, random rotation happens and new coordinates are saved in designated variables.
var point = [];
var pointlength = [];
var newXr = [];
var newYr = [];
var newZr = [];
//fragment loop
for (var m=1; m <= fragCount; m++) {
//random axis for fragment
var axisB = new Vector({ x: (Math.random())*2-1, y: (Math.random())*2-1, z: (Math.random())*2-1, w: 0 });
//fixed angle of 90 degrees (PI/2)
var angleC = Math.PI/2;
axisB.normalize();
var lenght_axisB = axisB.length();
//atom loop
for (var i=0; i < texts.length; i++) {
//check if atom is part of the fragment
if (Number(coorArrF[i]) == m) {
point[i] = new Vector({ x: newX[i]-centarx[m], y: newY[i]-centary[m], z: newZ[i]-centarz[m], w: 0 });
pointlength[i]= point[i].length();
point[i].normalize();
point[i].rotate(axisB, angleC);
point[i].scale(pointlength[i]);
newX[i]= point[i].x+centarx[m];
newY[i]= point[i].y+centary[m];
newZ[i]= point[i].z+centarz[m];
//finsih check
}
//finish atom loop
}
//finish fragment loop
}
After completing the new geometry, loop dynamics and all check are done.
It is interesting to note, that this is not a complete stochastic method, due to the fact that whole procedure starts from the relative positions of molecular fragments in input.
"Remembering" the initial positions severely limits the claim of totally random search of the chemical potential energy surface for a complete set of isomers. One can argue, that molecule simply wiggles from initial position. Hence the name of procedure.
KICK version
As stated above, original Saunders procedure randomly "kicks" atoms from their initial positions.To make a really relative initial position independent stochastic search, we have to choose randomly a new relative place from which "kick" starts should be chosen.
Using a bit of chemical logic here, we know that final result of the kick procedure should produce a sphere of fragments around the non-frozen fragment, which can be obtained from WIGGLE SCRIPT only if all fragments are put in the same spot and large R- distance parameter is supplied.
There are two approaches to do this:
1) VERSION find centers of all fragments and use origin as the relative positions for kick.
This is how WIGGLE SCRIPT could produce this result, but instead of manual input modification, it is done automatically. To do this, frozen fragment has to get it's center of fragment, and all coordinates have to recalculated relative to the center of the fragment. Essentially:
var relA = [];
var relX = [];
var relY = [];
var relZ = [];
//fragment loop
for (var m=0; m <= fragCount; m++) {
//atom loop
for (var i=0; i < texts.length; i++) {
//atom in fragment check
if (coorArrF[i] == m) {
var ada = coorA[i];
var xdx = (coorArrX[i]) - centarx[m] + 0.00001;
var ydy = (coorArrY[i]) - centary[m] + 0.00001;
var zdz = (coorArrZ[i]) - centarz[m] + 0.00001;
relA.push(ada);
relX.push(xdx);
relY.push(ydy);
relZ.push(zdz);
//atom in fragment check end
}
//atom loop end
}
//fragment loop end
}
Note: generation of
fragCountatomsInFrag[]centarX[]centarY[]centarZ[]Generation of random vector for each fragment.
var fragDx = [];
var fragDy = [];
var fragDz = [];
for (var m=1; m <= fragCount; m++) {
var vx = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vy = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vz = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vx2 = Math.pow(vx,2);
var vy2 = Math.pow(vy,2);
var vz2 = Math.pow(vz,2);
var vsum = vx2 + vy2 + vz2;
var vecnorm = Math.sqrt(vsum);
var vecsize = Math.random()*(param) + minDist;
var vectorX = (vx / vecnorm) * vecsize ;
var vectorY = (vy / vecnorm) * vecsize ;
var vectorZ = (vz / vecnorm) * vecsize ;
//VECTOR IN ARRAY OF VECTORS
var dx = vectorX;
var dy = vectorY;
var dz = vectorZ;
fragDx[m] = dx;
fragDy[m] = dy;
fragDz[m] = dz;
}
Now, to get a new coordinate, just add randomly generated vector to the relative atom positions.
var newA = [];
var newX = [];
var newY = [];
var newZ = [];
var newF = [];
for (var m=1; m <= fragCount; m++) {
for (var i=0; i < texts.length; i++) {
if (coorArrF[i] == m) {
var ada = coorA[i];
var xdx = relX[i] + fragDx[m];
var ydy = relY[i] + fragDy[m];
var zdz = relZ[i] + fragDz[m];
var fdf = coorArrF[i];
newA.push(ada);
newX.push(xdx);
newY.push(ydy);
newZ.push(zdz);
newF.push(fdf);
}
}
}
To finish everything, random rotation needs to be done on all non-frozen fragments, followed by distance calculation.
Although this is a acceptable approach, experience has shown that large R- distance parameter is needed to get reasonable structures. Better approach is the second version.
2) VERSION finds a random atom and use this as a starting point of the randomly generated vector.
This is achieved by:
var AllAtomsFragm = [];
AllAtomsFragm[0] = 0;
var RandomAtom = [];
RandomAtom[0] = 0;
for (var m=1; m <= fragCount; m++) {
AllAtomsFragm[m] = atomsInFrag[m-1] + AllAtomsFragm[m-1];
RandomAtom[m] = Math.floor(Math.random() * AllAtomsFragm[m]);
}
Note: generation of
fragCountatomsInFrag[]Generation of random vector.
var fragDx = [];
var fragDy = [];
var fragDz = [];
for (var m=1; m <= fragCount; m++) {
var vx = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vy = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vz = (Math.round(Math.random()) * 2 - 1)*Math.random()+0.0001;
var vx2 = Math.pow(vx,2);
var vy2 = Math.pow(vy,2);
var vz2 = Math.pow(vz,2);
var vsum = vx2 + vy2 + vz2;
var vecnorm = Math.sqrt(vsum);
var vecsize = Math.random()*(param) + minDist;
var vectorX = (vx / vecnorm) * vecsize ;
var vectorY = (vy / vecnorm) * vecsize ;
var vectorZ = (vz / vecnorm) * vecsize ;
//VECTOR IN ARRAY OF VECTORS
var dx = vectorX;
var dy = vectorY;
var dz = vectorZ;
fragDx[m] = dx;
fragDy[m] = dy;
fragDz[m] = dz;
ranX[m] = relX[RandomAtom[m]]+fragDx[m];
ranY[m] = relY[RandomAtom[m]]+fragDy[m];
ranZ[m] = relZ[RandomAtom[m]]+fragDz[m];
}
ranXranYranZNew coordinates are obtained as sum of relative atom coordinates with random point coordinates.
var ranX = [];
var ranY = [];
var ranZ = [];
var newX = [];
var newY = [];
var newZ = [];
var newA = [];
var newF = [];
for (var m=1; m <= fragCount; m++) {
ranX[m] = newX[RandomAtom[m]]+fragDx[m];
ranY[m] = newY[RandomAtom[m]]+fragDy[m];
ranZ[m] = newZ[RandomAtom[m]]+fragDz[m];
for (var i=0; i < texts.length; i++) {
if (coorArrF[i] == m) {
var ada = coorA[i];
var xdx = relX[i] + ranX[m];
var ydy = relY[i] + ranY[m];
var zdz = relZ[i] + ranZ[m];
var fdf = coorArrF[i];
newA.push(ada);
newX.push(xdx);
newY.push(ydy);
newZ.push(zdz);
newF.push(fdf);
}
}
}
To finish everything, random rotation needs to be done on all non-frozen fragments, followed by distance calculation.
This procedure can be retrofitted to be used for all atoms (as in original Saunders protocol), or in frozen fragment approach (fifth column of input is "0").
This is achieved in making one line in input reading:
var coorA = [];
var coorArrX = [];
var coorArrY = [];
var coorArrZ = [];
var coorArrF = [];
for (var i=0; i < texts.length; i++) {
var coor = texts[i].replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").split(" ");
//if not present as O, fifth column is set as a seperate entity and each atom is independently kicked
coor.push(i);
coorA[i] = coor[0];
coorArrX[i] = Number(coor[1]);
coorArrY[i] = Number(coor[2]);
coorArrZ[i] = Number(coor[3]);
coorArrF[i] = Number(coor[4]);
}
To conclude, we now have 5 distinct protocols to generate new random structures:
1) Original Saunders KICK version
2) Coalescence KICK version
3) Frozen fragment KICK version
4) WIGGLE version
5) KICK version
They all have strengths and weakness.
In our group we have over the years published several papers with different versions and procedures, including:
Original Saunders stochastic search procedure:
"Stochastic Search for Isomers of the sec-Butyl Cation
V. Vrček, O. Kronja, M. Saunders, J. Chem. Theory Comput. 2007, 3(3), 1223-1230.
Frozen fragment Saunders stochastic search procedure:
"Prereactive complexes in chlorination of benzene, triazine, and tetrazine: A quantum chemical study"
D. Šakić, V. Vrček, J. Phys. Chem. A 2012, 116(4), 1298-1306.
WIGGLE SCRIPT:
"Racemization of oxazepam and chiral 1,4-benzodiazepines. DFT study of the reaction mechanism in aqueous solution"
L. Hok, L. Božičević, H. Sremec, D. Šakić, V. Vrček, Org. Biomol. Chem. 2019, 17, 1471-1479.
KICK SCRIPT:
"A computational study of the chlorination and hydroxylation of amines by hypochlorous acid"
D. Šakić, M. Hanževački, D. M. Smith, V. Vrček, Org. Biomol. Chem. 2015, 13, 11740 - 11752.
I sincerely hope they will be useful to anyone that uses them!
Don't forget to cite and give credit!
Good luck!
In previous sections outline of the most vital javascript code was given. But for a script to be useful,
some html code was generated to wrap javascript and display results. First is the parameters form:
With following code:
<!--FORM GENERATION--> <form name="parameters"> Distance parameter: <input type="text" name="param" value="2.0"> (Angstrom) <br> Number of new files: <input type="text" name="counter" value="10"> <br> Minimal distance: <input type="text" name="minDist" value="0.1"> (Angstrom) <br> </form> <br> <!--FORM COLLECTION-->
Next, a text area for coordinates of structures can be inserted:
With code:
<textarea id="textinput" name="upload" cols="40" rows="10"> O -0.167787 1.645761 0.108747 0 H 0.613411 1.102620 0.113724 0 H -0.093821 2.209720 -0.643619 0 O 1.517569 -0.667424 -0.080674 1 H 1.989645 -1.098799 0.612047 1 H 0.668397 -1.091798 -0.139744 1 O 3.517569 -2.667424 -2.080674 2 H 3.989645 -3.098799 -1.588047 2 H 2.668397 -3.091798 -2.139744 2 </textarea>
And a button to set everything in motion:
Code:
Here function<button type="submit" class="generate" onClick="check()">GENERATE</button>
check() is activated by clicking submit. It is the major function in whichall prevously said procedures and loops are placed. In this example, important thing is added: a reporting
variable in which input when read is depostited. A lengthy code goes like this:
<script> function check(){ var lines = $('#textinput').val().split(/\n/); var texts = []; for (var i=0; i < lines.length; i++) { if (/\S/.test(lines[i])) { texts.push($.trim(lines[i])); } } var demo = []; demo += ("<br>"+"INPUT: "+"<br>"); for (var i=0; i < texts.length; i++) { var coor = texts[i].replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").replace(" "," ").split(" "); coor.push(0); coorA[i] = coor[0]; coorArrX[i] = Number(coor[1]); coorArrY[i] = Number(coor[2]); coorArrZ[i] = Number(coor[3]); coorArrF[i] = Number(coor[4]); demo += (coorA[i]+" "+coorArrX[i]+" "+coorArrY[i]+" "+coorArrZ[i]+" Frag: "+coorArrF[i]+"<br>"); } document.getElementById('demo').innerHTML = demo ; //this is the interface between javascript and html } </script>
Let's try this by clicking on this button:
Useful code:
//REPLACE FUNCTION html-break(br) with /n to convert reporting variable to downloadable version
function br2nl (str, replaceMode) {
var replaceStr = (replaceMode) ? "\n" : '';
return str.replace(/<\s*\/?br\s*[\/]?</gi, replaceStr);
}
//DOWNLOAD FUNCTION
function download(filename, text) {
var element = document.createElement('a');
element.setAttribute('href', 'data:text/plain;charset=utf-8,' + encodeURIComponent(text));
element.setAttribute('download', filename);
element.style.display = 'none';
document.body.appendChild(element);
element.click();
document.body.removeChild(element);
}
//VISUALIZATION (3Dmol.js; https://3dmol.csb.pitt.edu/index.html)
//*3Dmol.js: molecular visualization with WebGL
//N. Rego, D. Koes, Bioinformatics 2015, 31(8), 1322 - 1324.
//javascript part
var atoms = [{elem: 'O', x:-0.167787, y:1.645761, z:0.108747, bonds:[1,2]}, {elem: 'H', x:0.613411, y:1.102620, z:0.113724, bonds:[0]}, {elem: 'H', x:-0.093821, y:2.209720, z:-0.643619, bonds:[0]}, {elem: 'O', x:1.517569, y:-0.667424, z:-0.080674, bonds:[4,5]}, {elem: 'H', x:1.989645, y:-1.098799, z:0.612047, bonds:[3]}, {elem: 'H', x:0.668397, y:-1.091798, z:-0.139744, bonds:[3]},{elem: 'O', x:3.517569, y:-2.667424, z:-2.080674, bonds:[7,8]}, {elem: 'H', x:3.989645, y:-3.098799, z:-1.588047, bonds:[6]}, {elem: 'H', x:2.668397, y:-3.091798, z:-2.139744, bonds:[6]}];
$(function() {
let element = $('#container-01');
let config = { backgroundColor: 'white' };
let viewer = $3Dmol.createViewer( element, config );
viewer.setBackgroundColor(0xffffffff);
var m = viewer.addModel();
m.addAtoms(atoms);
m.setStyle({},{stick:{}});
viewer.zoomTo();
viewer.render();
});
<!--VIEWER 3Dmol.js--> <div id="container-02" class="mol-container1"></div> <style> .mol-container1 { width: 600px; height: 400px; position: relative; }</style>