All Downloads are FREE. Search and download functionalities are using the official Maven repository.

com.actelion.research.chem.phesaflex.EvaluableFlexibleOverlap Maven / Gradle / Ivy

There is a newer version: 2024.11.2
Show newest version
package com.actelion.research.chem.phesaflex;


import java.util.Map;

import com.actelion.research.chem.Coordinates;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.chem.alignment3d.PheSAAlignmentOptimizer.PheSASetting;
import com.actelion.research.chem.alignment3d.transformation.ExponentialMap;
import com.actelion.research.chem.alignment3d.transformation.Quaternion;
import com.actelion.research.chem.alignment3d.transformation.Rotation;
import com.actelion.research.chem.alignment3d.transformation.RotationDerivatives;
import com.actelion.research.chem.alignment3d.transformation.TransformationSequence;
import com.actelion.research.chem.alignment3d.transformation.Translation;
import com.actelion.research.chem.conf.BondRotationHelper;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.conf.TorsionDB;
import com.actelion.research.chem.forcefield.mmff.ForceFieldMMFF94;
import com.actelion.research.chem.optimization.Evaluable;
import com.actelion.research.chem.phesa.AtomicGaussian;
import com.actelion.research.chem.phesa.Gaussian3D;
import com.actelion.research.chem.phesa.MolecularVolume;
import com.actelion.research.chem.phesa.QuickMathCalculator;
import com.actelion.research.chem.phesa.ShapeVolume;
import com.actelion.research.chem.phesa.VolumeGaussian;
import com.actelion.research.chem.phesa.pharmacophore.pp.PPGaussian;
import com.actelion.research.chem.phesa.PheSAAlignment;



/**
 * @author JW, Oktober 2019
 * functionality for optimizing PheSA overlap (Pharmacophore+Shape) allowing for molecular flexibility
 */


public class EvaluableFlexibleOverlap implements Evaluable  {

	//private static final double SCALE = -250;
	//private static final double DELTA = -0.01;
	private static final double LAMBDA = 0.0625;
	private double e0 = 0.0;
	private Conformer fitConf;
	private StereoMolecule refMol;
	private PheSAAlignment shapeAlign;
	private boolean[] isHydrogen;
	private double[] v; // internal coordinates of the atoms
	private double[][] precalcPow;
	private double[] precalcExp;
    private double oAA;
    private double oAApp;
    private ForceFieldMMFF94 ff;
    private Map ffOptions;
    private PheSASetting settings;
	private Coordinates[] origCoords;
	private Coordinates[] cachedCoords; //used for gradient calculation: ligand coordinates with adjusted dihedral angles, but before rotation and translation
	private double[][] dRdvi1;
	private double[][] dRdvi2;
	private double[][] dRdvi3;
	private Coordinates origCOM;
	private BondRotationHelper torsionHelper;
    
	public EvaluableFlexibleOverlap(PheSAAlignment shapeAlign, StereoMolecule refMol, Conformer fitConf, BondRotationHelper torsionHelper,
			PheSASetting settings, boolean[] isHydrogen, Map ffOptions) {
		ForceFieldMMFF94.initialize(ForceFieldMMFF94.MMFF94SPLUS);
		this.ffOptions = ffOptions;
		this.shapeAlign = shapeAlign;
		this.fitConf = fitConf;
		this.isHydrogen = isHydrogen;
		this.settings = settings;
		this.refMol = refMol;
		this.torsionHelper = torsionHelper;
		init();
	}
	
	private void init() {
		ff = new ForceFieldMMFF94(fitConf.getMolecule(), ForceFieldMMFF94.MMFF94SPLUS, this.ffOptions);
		setInitialState();
		origCoords = new Coordinates[fitConf.getMolecule().getAllAtoms()];
		cachedCoords = new Coordinates[fitConf.getMolecule().getAllAtoms()];
		origCOM = new Coordinates();
		for(int a=0;a5) { //torsions
				if(v[i]>Math.PI) {
					v[i] -= 2*Math.PI;
				}
			}
			this.v[i] = v[i];
	
		}
		updateLigandCoordinates();
		shapeAlign.getMolGauss().update(fitConf);
	}
	
	public double[] getState(double[] v){
		for(int i=0;i0.0) {
						totalOverlap += atomOverlap;
						gradientPrefactor = atomOverlap*-2*refAt.getWidth()*fitAt.getWidth()/(refAt.getWidth()+fitAt.getWidth());
					}

				}
				grad[3*a] += (xj-xi)*gradientPrefactor;
				grad[3*a+1] += (yj-yi)*gradientPrefactor;
				grad[3*a+2] += (zj-zi)*gradientPrefactor;
				}

		
		}
		if(refMolGauss instanceof MolecularVolume) {
			for(VolumeGaussian refVG:((MolecularVolume)refMolGauss).getVolumeGaussians()){
			double xi = refVG.getCenter().x;
			double yi = refVG.getCenter().y;
			double zi = refVG.getCenter().z;
			for(AtomicGaussian fitAt:molGauss.getAtomicGaussians()){
				int a = fitAt.getAtomId();
				double atomOverlap = 0.0;
				double xj = v[3*a];
				double yj = v[3*a+1];
				double zj = v[3*a+2];
				double dx = xi-xj;
				double dy = yi-yj;
				double dz = zi-zj;
				double Rij2 = dx*dx + dy*dy + dz*dz;
				double alphaSum = refVG.getWidth() + fitAt.getWidth();
				double gradientPrefactor=0.0;
				if(Rij20.0) {
						totalOverlap += atomOverlap;
						gradientPrefactor = atomOverlap*-2*refVG.getWidth()*fitAt.getWidth()/(refVG.getWidth()+fitAt.getWidth());
					}

				}
				grad[3*a] += (xj-xi)*gradientPrefactor;
				grad[3*a+1] += (yj-yi)*gradientPrefactor;
				grad[3*a+2] += (zj-zi)*gradientPrefactor;
				}

		
		}
		}

		return totalOverlap; 
	
	}
	
	public double getFGValuePP() {

		ShapeVolume molGauss = shapeAlign.getMolGauss();

		ShapeVolume refMolGauss = shapeAlign.getRefMolGauss();

		/**
		 * derivative of ShapeOverlap with respect to the four elements of the quaternion and three elements of translation
		 * 
		 */ 
	    double totalOverlap = 0.0;
		for(PPGaussian refPP:refMolGauss.getPPGaussians()){
			double xi = refPP.getCenter().x;
			double yi = refPP.getCenter().y;
			double zi = refPP.getCenter().z;
			for(PPGaussian fitPP:molGauss.getPPGaussians()){
				double atomOverlap = 0.0;
				Coordinates fitCenterModCoord = fitPP.getCenter();
				double xj = fitCenterModCoord.x;
				double yj = fitCenterModCoord.y;
				double zj = fitCenterModCoord.z;
				double dx = xi-xj;
				double dy = yi-yj;
				double dz = zi-zj;
				double Rij2 = dx*dx + dy*dy + dz*dz;
				double alphaSum = refPP.getWidth() + fitPP.getWidth();
				if(Rij20.0) {
						double sim = refPP.getInteractionSimilarity(fitPP);
						atomOverlap *= sim;
						totalOverlap += atomOverlap;
					}

				}

				}
		
		}

		return totalOverlap; 
	
	}
	
	public double getFGValueShapeSelf(double[] grad, ShapeVolume molGauss,boolean rigid) {
		
		for(int i=0;i0.0) {
				gradientPrefactor = atomOverlap*-2*refAt.getWidth()*fitAt.getWidth()/(refAt.getWidth()+fitAt.getWidth());
			}

		}
			grad[3*b] += (2*xj-2*xi)*gradientPrefactor;
			grad[3*b+1] += (2*yj-2*yi)*gradientPrefactor;
			grad[3*b+2] += (2*zj-2*zi)*gradientPrefactor;
			
			return atomOverlap;
		}
	

	
	
	public double getFGValueSelfPP(ShapeVolume molVol,boolean rigid) {
		double xi,yi,zi,xj,yj,zj;

		/**
		 * derivative of ShapeOverlap with respect to Cartesian coordinates
		 */ 
	    double totalOverlap = 0.0;
	    for(PPGaussian refPP:molVol.getPPGaussians()){
			if(rigid) {
				xi = refPP.getCenter().x;
				yi = refPP.getCenter().y;
				zi = refPP.getCenter().z;
			}
			else {
				xi = refPP.getCenter().x;
				yi = refPP.getCenter().y;
				zi = refPP.getCenter().z;
			}
			for(PPGaussian fitPP:molVol.getPPGaussians()){
				double atomOverlap = 0.0;

				if(rigid) {
					xj = fitPP.getCenter().x;
					yj = fitPP.getCenter().y;
					zj = fitPP.getCenter().z;
				}
				else {
					xj = fitPP.getCenter().x;
					yj = fitPP.getCenter().y;
					zj = fitPP.getCenter().z;
				}
				double dx = xi-xj;
				double dy = yi-yj;
				double dz = zi-zj;
				double Rij2 = dx*dx + dy*dy + dz*dz;
				double alphaSum = fitPP.getWidth() + fitPP.getWidth();
				if(Rij20.0) {
						double sim = refPP.getInteractionSimilarity(fitPP);
						atomOverlap *= sim;
						totalOverlap += atomOverlap;
					}
				}
				}
		}
		return totalOverlap; 
	}
	
	public double[] getCartState(){
		double[] cartState = new double[3*fitConf.getMolecule().getAllAtoms()];
		for(int a=0;a




© 2015 - 2024 Weber Informatics LLC | Privacy Policy