com.actelion.research.chem.docking.LigandPose Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of openchemlib Show documentation
Open Source Chemistry Library
There is a newer version: 2024.12.1
package com.actelion.research.chem.docking;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Random;
import java.util.Set;
import java.util.stream.IntStream;

import com.actelion.research.calc.Matrix;
import com.actelion.research.chem.Coordinates;
import com.actelion.research.chem.StereoMolecule;
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.conf.TorsionDetail;
import com.actelion.research.chem.docking.InteractionTerm;
import com.actelion.research.chem.conf.torsionstrain.StatisticalTorsionPotential;
import com.actelion.research.chem.conf.torsionstrain.StatisticalTorsionTerm;
import com.actelion.research.chem.interactionstatistics.InteractionAtomTypeCalculator;
import com.actelion.research.chem.interactionstatistics.InteractionDistanceStatistics;
import com.actelion.research.chem.interactionstatistics.SplineFunction;
import com.actelion.research.chem.io.pdb.converter.MoleculeGrid;
import com.actelion.research.chem.phesa.Evaluable;
import com.actelion.research.chem.phesa.PheSAAlignment;
import com.actelion.research.chem.phesa.Quaternion;
import com.actelion.research.chem.potentialenergy.AngleConstraint;
import com.actelion.research.chem.potentialenergy.BondConstraint;
import com.actelion.research.chem.potentialenergy.PotentialEnergyTerm;
import com.actelion.research.chem.potentialenergy.TorsionConstraint;

import java.util.AbstractMap.SimpleEntry;

public class LigandPose implements Evaluable{
	
	private static double MOVE_AMPLITUDE = 2.0;
	private double BUMP_PENALTY = 500;
	private int BUMP_RADIUS = 3;
	private BondRotationHelper torsionHelper;
	private List ligAtomPairs; //separated by more than 3 bonds for internal strain
	//private double[] torsionValues;
	//private Quaternion rotation; //quaternion
	//private double[] translation;
	private double[] state; //coordinates
	private Conformer ligConf;
	private Conformer recConf;
	private StereoMolecule mol;
	private List ligStrain;
	private List interactionEnergy;
	private int[] ligAtomTypes;
	private int[] receptorAtomTypes;
	private Set receptorAtoms;
	private Random random;
	private MoleculeGrid grid;
	
	public LigandPose(Conformer ligConf, Conformer recConf, Set receptorAtoms, int[] receptorAtomTypes, 
			MoleculeGrid grid) {
		this.receptorAtomTypes = receptorAtomTypes;
		this.ligConf = ligConf;
		this.recConf = recConf;
		random = new Random(12345L);
		this.receptorAtoms = receptorAtoms;
		this.grid = grid;
		init();
	}
	
	private void init() {
		InteractionDistanceStatistics.getInstance().initialize();
		StatisticalTorsionPotential.getInstance().initialize();
		mol = ligConf.getMolecule();
		state = new double[3*mol.getAllAtoms()];
		ligStrain = new ArrayList();
		List ligAtomTypesList = new ArrayList<>();
		for(int a=0;a ligAtomTypes[e] = ligAtomTypesList.get(e));
		updateState();
		torsionHelper = new BondRotationHelper(mol);
		ligAtomPairs = new ArrayList();
		createStrainFunction();
		initiateInteractionTerms();
		
		
	}
	
	private double getBumpTerm() {
		double bumpTerm = 0.0;
		int[] gridSize = grid.getGridSize();
		for(Coordinates c : ligConf.getCoordinates()) {
			int[] gridC = grid.getGridCoordinates(c);
			int x = gridC[0];
			int y = gridC[1];
			int z = gridC[2];	
			if(x<0 || x>(gridSize[0]-3)) {
				bumpTerm = BUMP_PENALTY;
				break;
			}
			else if(y<0 || y>(gridSize[1])-3) {
				bumpTerm = BUMP_PENALTY;
				break;
			}
			else if(z<0 || z>(gridSize[1])-3) {
				bumpTerm = BUMP_PENALTY;
				break;	
				}
			}
		return bumpTerm;
		}
		

	
	private void createStrainFunction() {
		findLigAtomPairs();
		for(int[] pair : ligAtomPairs) {
			int at1 = pair[0];
			int at2 = pair[1];
			InteractionTerm term = InteractionTerm.create(ligConf,ligConf,at1,at2,ligAtomTypes,ligAtomTypes);
			if(term==null)
				continue;
			ligStrain.add(term);
		}
		//add torsions around rotatable bonds to term list and restrain the others to the current value, to prevent distortions
		for(int b=0;b 1) {
                for (int i=0; i> invalidPairs  = new HashSet>();
		StereoMolecule mol = ligConf.getMolecule();
		for(int a=0;a entry = a(a,aa) : new SimpleEntry<>(a,aa);
				invalidPairs.add(entry);
				for(int j=0;j(a,aaa) : new SimpleEntry<>(a,aaa);
					invalidPairs.add(entry);
					for(int k=0;k(a,aaaa) : new SimpleEntry<>(a,aaaa);
						invalidPairs.add(entry);
					}
				}
			}
		}
		for(int i=0;i entry =new SimpleEntry<>(i,j);
				if(invalidPairs.contains(entry))
					continue;
				else 
					ligAtomPairs.add(new int[] {i,j});
					
			}
		}
	}
		
		
	//constrain bonds that are not rotatable, constrain bond lengths and angles
	public double getFGValue(double[] gradient) {
		double energy = 0.0;
		for(PotentialEnergyTerm term : ligStrain) {
			energy+=term.getFGValue(gradient);
		}
		for(PotentialEnergyTerm term : interactionEnergy) {
			energy+=term.getFGValue(gradient);
		}
		energy+=getBumpTerm();
		return energy;
			
	}
	
	private void initiateInteractionTerms() {

		interactionEnergy = new ArrayList();

		for(int p : receptorAtoms) {
			for(int l=0;l0.0001) {
				Coordinates axis = rot.scale(1.0/angle);
				q = new Quaternion(axis,angle);
			}
			Matrix m = q.getRotMatrix();
			Coordinates com = DockingUtils.getCOM(ligConf);
			ligConf.translate(-com.x, -com.y, -com.z);
			PheSAAlignment.rotateMol(ligConf, m);
			ligConf.translate(com.x, com.y, com.z);
			
		}
		else if(num==2) {
			double targetTorsion = 360*random.nextDouble();
			int bond = random.nextInt(torsionHelper.getRotatableBonds().length);
			//calculate actual torsion, get difference to desired torsion and take
			//code snippet from torsionHelper to rotate the atoms
			int[] torsionAtoms = torsionHelper.getTorsionAtoms()[bond];
			Coordinates c1 = ligConf.getCoordinates(torsionAtoms[0]);
			Coordinates c2 = ligConf.getCoordinates(torsionAtoms[1]);
			Coordinates c3 = ligConf.getCoordinates(torsionAtoms[2]);
			Coordinates c4 = ligConf.getCoordinates(torsionAtoms[3]);
			double currentTorsion = Coordinates.getDihedral(c1, c2, c3, c4);
			if(currentTorsion<0)
				currentTorsion+=2*Math.PI;
			currentTorsion = 360.0*currentTorsion/2*Math.PI;
			double rotateBy = targetTorsion-currentTorsion;
			torsionHelper.rotateSmallerSide(bond, rotateBy);
			
		}
		updateState();
		
	}
		
}