#! /bin/awk -f
#
#	Generate a superhelix
#
#	This program reads one line from standard input containing 10 numbers (Crick parameters) which have default values:
#	R0 = 4.809		# superhelical radius (Angstroms)
#	pitch = -109.642	# superhelical pitch (residues/superturn)
#       phi = 13.214            # phase of helix vs superhelix (degrees)
#	helicies = 2		# number of helicies in oligomer
#	nres = 33		# number of aa residues in each helix
#       R1 = 2.269              # intra-helical radius (Angstroms)
#       d  = 1.497              # intra-helical rise/residue  (Angstroms)
#       w1 = 101.808            # intra-helical twist/residue (degrees)
#       W = 0                   # rotation of superhelix (degrees)
#	zc = 0			# z-offset from origin (Angstroms)
#
#
#	You may list less than 10 numbers and/or replace one of these numbers with a "-" to use the default value
#
#	Example of how to run it:
#	echo "4.8 -110 - 3 33 - - - 10 -20 " | awk -f supertwist.awk > test.pdb
#
#
BEGIN{
}

{
    R0    = $1;
    pitch = $2;
    phi   = $3;
    helicies = $4;
    nres  = $5;
    R1    = $6;
    d     = $7;
    w1    = $8;
    W     = $9;
    zc    = $10;

    R2 = $11
    phi2 = $12
    z2   = $13

    supertwist(R0,pitch,phi,helicies,nres,R1,d,w1,zc,W)
}


function supertwist(R0,pitch,phi,helicies,nres,R1,d,w1,zc,W)
{
    PI = 2*atan2(1,0)

    # "-" indicates default value (GCN4)
    if(R0 == "-") R0 = ""
    if(pitch == "-") pitch = ""
    if(helicies == "-") helicies = ""
    if(nres == "-") nres = ""
    if(phi=="-") phi = ""
    if(zc == "-") zc = ""
    if(W == "-") W = ""
    if(R1== "-") R1 = ""
    if(d == "-") d  = ""
    if(w1== "-") w1 = ""

    # sensible defaults? (classic gcn4 dimer)
    if(R0 == "")        R0 = 4.809	# superhelical radius (Angstroms)
    if(pitch == "") pitch = -109.642	# superhelical pitch (residues/superturn)

    if(helicies == "") helicies = 2	# number of helicies in oligomer
    if(nres == "")         nres = 33	# number of aa residues in each helix

    if(phi=="")           phi = 13.214	# phase of helix vs superhelix (degrees)
    if(zc == "")           zc = 0	# z-offset from origin (Angstroms)
    if(W == "") 	    W = 0	# rotation of superhelix (degrees)

    if(R1== "") R1 = 2.269		# helical radius (Angstroms)
    if(d == "") d  = 1.497		# helical rise/residue  (Angstroms)
    if(w1== "") w1 = 101.808		# helical twist/residue (degrees)

    # convert to radians
    phi   *= PI/180
    W     *= PI/180
    w1    *= PI/180

    # superhelical twist/residue
    if(pitch == "0") w0 = 0
    if(pitch+0 != 0) w0 = 2*PI/pitch

    # sanity check
    if(w0^2 > (d/R0)^2) 
    {
	# supertwist is tighter than helical twist
	print "REMARK WARNING: helical pitch changed to", (2*PI)/w0
	w0=d/R0
    }
    if(helicies < 1) helicies = 1
    if(nres < 0) nres = -nres
    if(nres < 1) nres = 1

    # helical crossing angle
    alpha = asin(R0*w0/d)

    for(i=0;i<helicies;++i) 
    {
	chain = sprintf("%c", i+65)
	for(t=1;t<=nres;++t)
	{
	    ++Atomnum
		
	    # determine position of N atom
	    RN   = 0.661*R1;
	    phiN = phi-0.265*w1;
	    dN   = -0.584*d
#RN= R2*R1;
#phiN= phi +phi2*w1
#dN= z2*d;
	    N["X"] = RN * cos(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phiN) \
              - RN * sin(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phiN) * cos(alpha) \
	      + R0 * cos(w0*t + W + 2*PI*i/helicies)

	    N["Y"] = RN * sin(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phiN) \
	      + RN * cos(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phiN) * cos(alpha) \
	      + R0 * sin(w0*t + W + 2*PI*i/helicies)

	    N["Z"] = -RN*sin(w1*t + phiN)*sin(alpha) + t*d*cos(alpha) + zc +dN


	    # determine position of C atom
	    RC   = 0.745*R1;
	    phiC = phi+0.272*w1;
	    dC   = 0.689*d;
#RC= R2*R1;
#phiC= phi +phi2*w1
#dC= z2*d;
	    C["X"] = RC * cos(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phiC) \
              - RC * sin(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phiC) * cos(alpha) \
	      + R0 * cos(w0*t + W + 2*PI*i/helicies)

	    C["Y"] = RC * sin(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phiC) \
	      + RC * cos(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phiC) * cos(alpha) \
	      + R0 * sin(w0*t + W + 2*PI*i/helicies)

	    C["Z"] = -RC*sin(w1*t + phiC)*sin(alpha) + t*d*cos(alpha) + zc +dC


	    # determine position of O atom
	    RO   = 0.897*R1;
	    phiO = phi+0.223*w1;
	    dO   = 1.447*d;
#RO= R2*R1;
#phiO= phi +phi2*w1
#dO= z2*d;
	    O["X"] = RO * cos(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phiO) \
              - RO * sin(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phiO) * cos(alpha) \
	      + R0 * cos(w0*t + W + 2*PI*i/helicies)

	    O["Y"] = RO * sin(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phiO) \
	      + RO * cos(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phiO) * cos(alpha) \
	      + R0 * sin(w0*t + W + 2*PI*i/helicies)

	    O["Z"] = -RO*sin(w1*t + phiO)*sin(alpha) + t*d*cos(alpha) + zc +dO


	    # determine position of CA atom
	    CA["X"] = R1 * cos(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phi) \
              - R1 * sin(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phi) * cos(alpha) \
	      + R0 * cos(w0*t + W + 2*PI*i/helicies)

	    CA["Y"] = R1 * sin(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phi) \
	      + R1 * cos(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phi) * cos(alpha) \
	      + R0 * sin(w0*t + W + 2*PI*i/helicies)

	    CA["Z"] = -R1*sin(w1*t + phi)*sin(alpha) + t*d*cos(alpha) + zc


	    # determine position of CB atom
	    RB   = 1.460*R1;
	    phiB = phi+0.172*w1;
	    dB   = -0.492*d;
#RB= R2*R1;
#phiB= phi +phi2*w1
#dB =  z2*d;
	    CB["X"] = RB * cos(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phiB) \
              - RB * sin(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phiB) * cos(alpha) \
	      + R0 * cos(w0*t + W + 2*PI*i/helicies)

	    CB["Y"] = RB * sin(w0*t + W + 2*PI*i/helicies) * cos(w1*t + phiB) \
	      + RB * cos(w0*t + W + 2*PI*i/helicies) * sin(w1*t + phiB) * cos(alpha) \
	      + R0 * sin(w0*t + W + 2*PI*i/helicies)

	    CB["Z"] = -RB*sin(w1*t + phiB)*sin(alpha) + t*d*cos(alpha) + zc +dB

	    if(! restyp) restyp = "ALA"
	    if(! occ)    occ = 1
	    if(! bfac)   Bfac = 80

	    printf "%s", sprint_atom(N,"N",restyp,t);
#	    occ = 0
	    printf "%s", sprint_atom(CA,"CA",restyp,t);
	    occ = 1
	    printf "%s", sprint_atom(C,"C",restyp,t);
	    printf "%s", sprint_atom(O,"O",restyp,t);
	    printf "%s", sprint_atom(CB,"CB",restyp,t);
	}
    }

}

################################################################################
#
#       sprint_atom(atom, name, restyp, resnum, occ, Bfac)
#
#       Function for creating a standard PDB line
#
################################################################################
function sprint_atom(atom, _name, _restyp, _resnum, _occ, _Bfac) {
# defaults (global)
if(! name)     name = "CA"
if(! restyp)   restyp = "ALA"
if(resnum=="") resnum = 1
if(! chain)    chain = ""
if(occ=="")    occ = 1
if(Bfac=="")   Bfac = 20

if(_name)   name   = toupper(_name)
if(_restyp) restyp = toupper(_restyp)
if(_resnum) resnum = toupper(_resnum)
if(_occ)    occ = _occ
if(_Bfac)   Bfac = _Bfac

if(resnum ~ /^[A-Z]/) chain = substr(resnum,1,1)
while(resnum ~ /^[A-Z]/) resnum = substr(resnum,2)
++atoms_printed

entry=sprintf("ATOM %6d  %-3s %3s %1s%4d     %7.3f %7.3f %7.3f %5.2f%6.2f\n",\
atoms_printed,name,restyp,chain,resnum,atom["X"],atom["Y"],atom["Z"],occ,Bfac);

return entry
}


################################################################################
#
#	next_atom(atom1,atom2,atom3, chi, angle, bond)
#
#  	Function for getting "new atom" xyz coordinates using:
#
#	three reference atoms (defining two "previous" bonds)
#	two   angles (the bond angle, and the chi torsion angle)
#	one   distance (length of the "new" bond)
#
#
#     O -atom1                     O - new_atom
#      \                          /
#       \                        /
#        \           angle      / - bond (1.54A)
#         \               \__  /
#          \   chi(= 0)   /   /
#    atom2- O -------------- O -atom3
#
# atom1, atom2, atom3, are 3-membered arrays ["X","Y","Z"]
# new atom returned in new_atom["X","Y","Z"]
#
# note: atom1 serves only to define chi=0, it can really be anywhere
################################################################################
#
#                                   |-> optional
function next_atom(atom1,atom2,atom3, chi, angle, bond) {


    # default bond angle and length of new bond
    if(chi   == "") chi   = 0
    if(angle == "") angle = 109.5
    if(bond == "")
    {
	bond = 1.54
	# assume some double-bond character in non-tetrahedral bonds
	if(angle != 109.5) bond = 1.4
    }
    

    # compute components of "new_atom"-"atom3" vector, relative to chi==0 defined by "bond1"
     axis_component["length"] = bond*sin(3.1415927*(angle-90)/180)
     chi0_component["length"] = -bond*cos(3.1415927*(angle-90)/180)*cos(3.1415927*(chi)/180)
    chi90_component["length"] = -bond*cos(3.1415927*(angle-90)/180)*sin(3.1415927*(chi)/180)

    # now we know the coordinates of the new atom vector in "local" coordinates
    # we need to construct a basis for converting them to "global" coordinates

    # vector subtration of atoms lying on rotation axis
    axis["X"] = atom3["X"]-atom2["X"]
    axis["Y"] = atom3["Y"]-atom2["Y"]
    axis["Z"] = atom3["Z"]-atom2["Z"]

    # vector subtraction of atoms defining "zero" rotation around the axis
    bond1["X"] = atom2["X"]-atom1["X"]
    bond1["Y"] = atom2["Y"]-atom1["Y"]
    bond1["Z"] = atom2["Z"]-atom1["Z"]

    # protect against singular vectors
    if(((axis["X"]^2 + axis["Y"]^2 + axis["Z"]^2) == 0)||((bond1["X"]^2 + bond1["Y"]^2 + bond1["Z"]^2)==0))
    {
	new_atom["X"] = 0;
	new_atom["Y"] = 0;
	new_atom["Z"] = 0; 
	return 0   
    }

    # normalize the "axis" vector
    axis["length"]  = sqrt( (axis["X"]*axis["X"])  + ( axis["Y"]*axis["Y"])  + (axis["Z"]*axis["Z"]));
    axis["X"] = axis["X"]/axis["length"]
    axis["Y"] = axis["Y"]/axis["length"]
    axis["Z"] = axis["Z"]/axis["length"]
    axis["length"] = 1

    # compute amount of "bond1" that needs to be removed in order to orthogonalize it to "axis"
    bond1_dot_axis  =     ( (axis["X"]*bond1["X"]) +  (axis["Y"]*bond1["Y"]) + (axis["Z"]*bond1["Z"]) ) # Dot product

    # subtract the projection from bond1
    bond1["X"] = bond1["X"] - bond1_dot_axis*axis["X"]
    bond1["Y"] = bond1["Y"] - bond1_dot_axis*axis["Y"]
    bond1["Z"] = bond1["Z"] - bond1_dot_axis*axis["Z"]

    # normalize the "bond1" vector to make the "chi0" reference vector
    bond1["length"] = sqrt((bond1["X"]*bond1["X"]) + (bond1["Y"]*bond1["Y"]) +(bond1["Z"]*bond1["Z"]) );
    chi0["X"] = bond1["X"]/bond1["length"]
    chi0["Y"] = bond1["Y"]/bond1["length"]
    chi0["Z"] = bond1["Z"]/bond1["length"]
    chi0["length"] = 1

    # now make the "other" basis vector to complete the right-handed, 
    # orthonormal basis, using a cross-product
    chi90["X"] = axis["Y"] * chi0["Z"] - axis["Z"] * chi0["Y"];
    chi90["Y"] = axis["Z"] * chi0["X"] - axis["X"] * chi0["Z"];
    chi90["Z"] = axis["X"] * chi0["Y"] - axis["Y"] * chi0["X"];
    
    # we now have three unit vectors forming a basis of rotation about the "atom2-atom3" bond

    # the "axis" component of the new atom, offset from atom3, will be 0.514A out along "axis"
    axis_component["X"]  = axis_component["length"]*axis["X"]
    axis_component["Y"]  = axis_component["length"]*axis["Y"]
    axis_component["Z"]  = axis_component["length"]*axis["Z"]

    # apply the "x-y" values from the given chi angle
    chi0_component["X"]  = chi0_component["length"]*chi0["X"]
    chi0_component["Y"]  = chi0_component["length"]*chi0["Y"]
    chi0_component["Z"]  = chi0_component["length"]*chi0["Z"]

    chi90_component["X"] = chi90_component["length"]*chi90["X"]
    chi90_component["Y"] = chi90_component["length"]*chi90["Y"]
    chi90_component["Z"] = chi90_component["length"]*chi90["Z"]

    
    # now generate a "new" atom, in the original coordinate system
    new_atom["X"] = atom3["X"] + axis_component["X"] + chi0_component["X"] + chi90_component["X"]
    new_atom["Y"] = atom3["Y"] + axis_component["Y"] + chi0_component["Y"] + chi90_component["Y"]
    new_atom["Z"] = atom3["Z"] + axis_component["Z"] + chi0_component["Z"] + chi90_component["Z"]

    return 1
}

function asin(_sin)
{
    #print "atan2("_sin","sqrt(1-_sin*_sin)")="atan2(_sin,sqrt(1-_sin*_sin))
    return atan2(_sin,sqrt(1-_sin*_sin))
}