c Samples ligand orientations around a covalent bond. c Given as input the two atom locations of the covalent bond c (the latter previuosly generated by cov_sampler), c an array of the ligand corrdinates, and the index in the array c for the ligand covalent modification site (ligi1) and it's c prior atom (ligi2) aligns the ligand in correct bond angle, c and dihedrals. c NOTE: * fi1 and fi2 are relative to prior ligand position. c * The a param is absolute if na is true and relative c if it is false - this is confusing behaviour, but at c least it's documented. c NL 01/2012 subroutine cov_lig (p1,p2,ligi1,ligi2,ligi3,ligi4, & lig,ligsize,a,fi1,fi2,sp,na) implicit none c (receptor) c * p1 c \_ c fi1 /_\/ ______ c \--- / \ c \ a \ /\ ligi2 rest \ c *-----|-----* of > c p2&ligi1 \/ \ lig / c (lig) fi2 \______/ c input variables real p1(3) !receptor's covalent attachment site real p2(3) !ligand covalent attachment site integer ligi1,ligi2,ligi3,ligi4 !indices of ligand atoms cov,n1-3 integer ligsize !number of ligand atoms real lig(3,ligsize) !ligand's atoms coordinates real a,fi1,fi2 !ligand's placement params see fig. integer sp !indicated the hybridization geometry of the cov atom logical na !if true 'a' is absolute, !if false 'a' is relative c local variables integer i,j !loop variables real l1(3), l2(3) !ligand atoms real cov(3), lig_bond(3), cross_cl(3) !local vectors real lig31(3),lig32(3),lig42(3) !local vectors real lig32x42(3), lig42x32(3) !local vectors real tr(3) !translation vector to take ligand to p2 real rot1(3,3) !rot mat. to align the cov bond with ligand real rot2(3,3) !rot mat. to create bond angle: a real rot3(3,3) !rot mat. to create dihedral: fi1 real rot4(3,3) !rot mat. to create dihedral: fi2 real rot_sp3(3,3) !rot mat. to align sp3 hybrid. real rot_sp2_1(3,3) !first rot. mat. for sp2 real rot_sp2_2(3,3) !second rot. mat. for sp2 c constants REAL, PARAMETER :: Pi = 3.1415927 REAL, PARAMETER :: deg_to_rad = 0.0174532925 !debug real a,b c main code c ========= c translate the ligand from its current arbitrary location c so that ligi1 coinside with p2 tr(1)=p2(1)-lig(1,ligi1); tr(2)=p2(2)-lig(2,ligi1); tr(3)=p2(3)-lig(3,ligi1); !apply the translation do j=1,ligsize lig(1,j)=lig(1,j)+tr(1); lig(2,j)=lig(2,j)+tr(2); lig(3,j)=lig(3,j)+tr(3); enddo !vector representing the covalent bond cov=p2-p1; !calculate the initial ligand orientation (to be applied for any new !covalent attachment point creates the first bond angle (by rotation !matrix to align the first ligand bond and cov bond at 0 deg !and then orients by sp2 or sp3 geometry if (na) then !ligand atoms coordinates l1(1)=lig(1,ligi1); l1(2)=lig(2,ligi1); l1(3)=lig(3,ligi1); l2(1)=lig(1,ligi2); l2(2)=lig(2,ligi2); l2(3)=lig(3,ligi2); lig_bond=l2-l1; !cross product of cov and lig bonds, to find the axis !around which to create the bond angle. cross_cl(1)=lig_bond(2)*cov(3)-lig_bond(3)*cov(2) cross_cl(2)=lig_bond(3)*cov(1)-lig_bond(1)*cov(3) cross_cl(3)=lig_bond(1)*cov(2)-lig_bond(2)*cov(1) call rotate(lig_bond,cov,cross_cl,rot1) !calculate rotation matrix to create bond angle 'a' call rotate_ang(Pi-(deg_to_rad*a),cross_cl,rot2) endif ! na check !calculate rotation matrix to create dihedral 'fi1' !NOTE IMPORTANT - this only increments the angle by fi1 degs call rotate_ang(deg_to_rad*fi1,cov,rot3) !translate ligand to local 'origin' do i=1,3 do j=1,ligsize lig(i,j)=lig(i,j)-p2(i) enddo enddo !apply rotations if (na) then lig=matmul(rot1,lig) lig=matmul(rot2,lig) !recalculate lig vector after initial rotations lig_bond(1)=lig(1,ligi2)-lig(1,ligi1); lig_bond(2)=lig(2,ligi2)-lig(2,ligi1); lig_bond(3)=lig(3,ligi2)-lig(3,ligi1); !for sp3 hyb. calculate the normal to the plane defined by the ! three neighbours and align it to the cov bond. if (sp .eq. 3) then lig32(1)=lig(1,ligi3)-lig(1,ligi2); lig32(2)=lig(2,ligi3)-lig(2,ligi2); lig32(3)=lig(3,ligi3)-lig(3,ligi2); lig42(1)=lig(1,ligi4)-lig(1,ligi2); lig42(2)=lig(2,ligi4)-lig(2,ligi2); lig42(3)=lig(3,ligi4)-lig(3,ligi2); !calculate both normals to the planes lig32x42(1)=lig32(2)*lig42(3)-lig32(3)*lig42(2) lig32x42(2)=lig32(3)*lig42(1)-lig32(1)*lig42(3) lig32x42(3)=lig32(1)*lig42(2)-lig32(2)*lig42(1) lig42x32=-1*lig32x42 !the normals to the plane are going through !the origin - i.e. ligi1, so, to select a normal, see which !is closer to e.g. ligi2... !check which is closer to ligi1 to decide which to align to cov if ( & (((lig(1,ligi2)-lig32x42(1))*(lig(1,ligi2)-lig32x42(1)))+ & ((lig(2,ligi2)-lig32x42(2))*(lig(2,ligi2)-lig32x42(2)))+ & ((lig(3,ligi2)-lig32x42(3))*(lig(3,ligi2)-lig32x42(3)))) < & (((lig(1,ligi2)-lig42x32(1))*(lig(1,ligi2)-lig42x32(1)))+ & ((lig(2,ligi2)-lig42x32(2))*(lig(2,ligi2)-lig42x32(2)))+ & ((lig(3,ligi2)-lig42x32(3))*(lig(3,ligi2)-lig42x32(3)))) & ) then !vector representing the orientation of the normal to the plane !should be aligned with the cov_bond vector !if (1 .ne. 0) then call rotate(lig32x42,cov,lig_bond,rot_sp3) else call rotate(lig42x32,cov,lig_bond,rot_sp3) endif lig=matmul(rot_sp3,lig) endif ! sp==3 !for sp2 hyb. if (sp .eq. 2) then !calculate the lig3 to lig1 vector lig31(1)=lig(1,ligi3)-lig(1,ligi1); lig31(2)=lig(2,ligi3)-lig(2,ligi1); lig31(3)=lig(3,ligi3)-lig(3,ligi1); !recalculate the normal to the plane created by cov bond and lig 21 call rotate(lig31,cross_cl,lig_bond,rot_sp2_1) call rotate_ang(deg_to_rad*90.0,lig_bond,rot_sp2_2) lig=matmul(rot_sp2_1,lig) lig=matmul(rot_sp2_2,lig) endif endif ! check na - first orientation !apply fi1 rotation lig=matmul(rot3,lig) !translate ligand back to p2 do i=1,3 do j=1,ligsize lig(i,j)=lig(i,j)+p2(i) enddo enddo return end