// put_tmp is a helper function to simplify the code somewhat -
// all it does is output the proper elements to a temporary
// pdb file.

int put_tmp( float x, float y, float z, int count, int rescount,
string atomname, string resname, file outfile, int nres ) {

// check if 5' end

if( rescount == 1 || rescount == ( nres / 2 ) + 1 ) 
	if( atomname =~ "P" )
			return 1;

if( atomname == "O1'" )
	atomname = "O4'";

// I have no idea why this is necessary.. would be just as easy
// to change O1' to O4' in all data files, but they must have had
// some reason for naming it O1' in the data.

fprintf(outfile, "ATOM  %5d %-4s %3s%6d    %8.3lf%8.3lf%8.3lf\n",
	count, atomname, resname, rescount, x, y, z );

return 1;

};

// fd_helix is a function designed to build helices of varying types
// from fiber diffraction data.  The current available types are:

// arna		Right Handed A-RNA (Arnott)
// aprna		Right Handed A-PRIME RNA (Arnott)
// lbdna		Right Handed B-DNA (Langridge)
// abdna		Right Handed B-DNA (Arnott)
// sbdna		Left  Handed B-DNA (Sasisekharan)
// adna		Right Handed A-DNA (Arnott)

// fd_helix will create a temporary file in your working directory
// called nab_tmp.pdb.  This file contains the helix _before_ the
// addition of hydrogens.  This is ordinarily deleted, but could be
// examined if things go wrong.


molecule fd_helix( string helix_type, string seq, string acid_type )
{

// ------------------------------------------------------------------ //

//	Data Segment

molecule m;

// These arrays are used to store the r, phi, zz values of the
// various atom types within each residue:

float ade_r[ hashed ];
float ade_phi[ hashed ];
float ade_zz[ hashed ];

float gua_r[ hashed ];
float gua_phi[ hashed ];
float gua_zz[ hashed ];

float thy_r[ hashed ];
float thy_phi[ hashed ];
float thy_zz[ hashed ];

float ura_r[ hashed ];
float ura_phi[ hashed ];
float ura_zz[ hashed ];

float cyt_r[ hashed ];
float cyt_phi[ hashed ];
float cyt_zz[ hashed ];

// These are simply temp variables used to aid in filling the arrays.

float temp_r, temp_phi, temp_zz;

// These are used to calculate and output coordinates to the temporary pdb file.

float x, y, z, yyr, xrad;

// height values are angstroms, rotation; values are degrees.

float current_height, current_rotation;
float height_increment, rotation_increment;

float hxht[ hashed ];
float hxrep[ hashed ];

string tempname, cseq, fullseq, buffer, restype;
string temp, amberhome;
string resout;

int nres, nresh, i, hxmul, count, chain, begin, end;

file infile, outfile;

// -----------------------------------------------------------

// 	Code Segment
hxht[ "arna" ] = 2.81; hxht[ "aprna" ] = 3.00; hxht[ "lbdna" ] = 3.38;
hxht[ "abdna" ] = 3.38; hxht[ "sbdna" ] = -3.38; hxht[ "adna" ] = 2.56;

hxrep["arna"] = 32.7; hxrep["aprna"] = 30.0; hxrep["lbdna"] = 36.0;
hxrep["abdna"] = 36.0; hxrep["sbdna"] = 36.0; hxrep["adna"] = 32.7;

temp = wc_complement( seq, acid_type + "amber94.rlb", acid_type );
cseq = "";
for( i = length( temp ); i >= 1; i-- ) {
	cseq += substr( temp, i, 1) ;
}
fullseq = seq + cseq;

nresh = length( seq );

nres = length( fullseq );

if( !( amberhome = getenv( "AMBERHOME" ) ) ){
    fprintf( stderr, "AMBERHOME not defined.\n" );
    exit( 1 );
}

temp = amberhome + "/dat/fd_data/" + helix_type + ".dat";
printf("Data file: %s\n", temp);

infile = fopen( temp, "r" );
if( infile == NULL ){
	fprintf( stderr, "Unable to open data file %s; exiting\n", temp );
	exit(1);
}

outfile = fopen( "nab_tmp.pdb", "w" );

// Read lines from data file and store r, phi, zz values as appropriate
// residue type.

while( buffer = getline( infile ) ) {
	sscanf( buffer, "%s %lf %lf %lf %s", tempname, temp_r, temp_phi,
		temp_zz, restype );

	if( restype == "A" || restype == "a" ) {
		ade_r[ tempname ] = temp_r;
		ade_phi[ tempname ] = temp_phi;
		ade_zz[ tempname ] = temp_zz;
	}
	else if( restype == "G" || restype == "g" ) {
		gua_r[ tempname ] = temp_r;
		gua_phi[ tempname ] = temp_phi;
		gua_zz[ tempname ] = temp_zz;
	}
	else if( restype == "T" || restype == "t" ) {
		thy_r[ tempname ] = temp_r;
		thy_phi[ tempname ] = temp_phi;
		thy_zz[ tempname ] = temp_zz;
	}
	else if( restype == "U" || restype == "u" ) {
		ura_r[ tempname ] = temp_r;
		ura_phi[ tempname ] = temp_phi;
		ura_zz[ tempname ] = temp_zz;
	}
	else if( restype == "C" || restype == "c" ) {
		cyt_r[ tempname ] = temp_r;
		cyt_phi[ tempname ] = temp_phi;
		cyt_zz[ tempname ] = temp_zz;
	}

}

height_increment = hxht[ helix_type ]; 
rotation_increment = hxrep[ helix_type ]; 
current_height = 0;
current_rotation = 0;
count = 0;

// 	Here we build the actual molecule - it is output to a temporary
//	pdb file in order to allow the addition of hydrogens as a final
// 	stage.

for( chain = 1; chain <= 2; chain++ ) {
	
	if(chain == 1) {
		begin = 1;
		end = nresh;
		hxmul = -1;
	}

	else if( chain == 2) {
		begin = nresh + 1;
		end = nres;
		hxmul = 1;
	}
	for( i = begin; i <= end; i++ ) {
		restype = substr( fullseq, i, 1 );
		if( restype == "A" || restype == "a" ) {
			resout = "DA"; if( acid_type == "rna" ) resout = "A";
			for( tempname in ade_r ) {
				count++;
				yyr = (hxmul * ade_phi[ tempname ] + current_rotation);
				xrad = ade_r[ tempname ];
				x = xrad * cos(yyr);
				y = xrad * sin(yyr);
				z = hxmul * ade_zz[ tempname ] + current_height;
				put_tmp( x, y, z, count, i, tempname, resout, outfile, nres);
			}	
		}

		else if( restype == "G" || restype == "g" ) {
			resout = "DG"; if( acid_type == "rna" ) resout = "G";
			for( tempname in gua_r ) {	
				count++;
				yyr = (hxmul * gua_phi[ tempname ] + current_rotation);
				xrad = gua_r[ tempname ];
				x = xrad * cos(yyr);
				y = xrad * sin(yyr);
				z = hxmul * gua_zz[ tempname ] + current_height;
				put_tmp( x, y, z, count, i, tempname, resout, outfile, nres);
			}	
		
		}

		else if( restype == "T" || restype == "t" ) {
			resout = "DT";
			for( tempname in thy_r ) {	
				count++;
				yyr = (hxmul * thy_phi[ tempname ] + current_rotation);
				xrad = thy_r[ tempname ];
				x = xrad * cos(yyr);
				y = xrad * sin(yyr);
				z = hxmul * thy_zz[ tempname ] + current_height;
				put_tmp( x, y, z, count, i, tempname, resout, outfile, nres);
			}	
		
		}

		else if( restype == "U" || restype == "u" ) {
			resout = "U";
			for( tempname in ura_r ) {	
				count++;
				yyr = (hxmul * ura_phi[ tempname ] + current_rotation);
				xrad = ura_r[ tempname ];
				x = xrad * cos(yyr);
				y = xrad * sin(yyr);
				z = hxmul * ura_zz[ tempname ] + current_height;
				put_tmp( x, y, z, count, i, tempname, resout, outfile, nres);
			}	
			
		}

		else if( restype == "C" || restype == "c" ) {
			resout = "DC"; if( acid_type == "rna" ) resout = "C"; 
			for( tempname in cyt_r ) {	
				count++;
				yyr = (hxmul * cyt_phi[ tempname ] + current_rotation);
				xrad = cyt_r[ tempname ];
				x = xrad * cos(yyr);
				y = xrad * sin(yyr);
				z = hxmul * cyt_zz[ tempname ] + current_height; 
				put_tmp( x, y, z, count, i, tempname, resout, outfile, nres);
			}	
		
		}
		
		// Increase unit twist and height

		current_height += height_increment;
		current_rotation += rotation_increment;
	}


height_increment = -height_increment;
rotation_increment = -rotation_increment;

current_rotation += rotation_increment;
current_height += height_increment;

if(chain == 1)
	fprintf( outfile, "TER\n" );	// Need TER card for getpdb_prm to
					// function properly

}

fclose( infile );

// Close outfile in order to flush output stream.
// Otherwise nab_tmp.pdb would be incomplete when
// read by getpdb_prm

fclose( outfile );

m = getpdb_prm( "nab_tmp.pdb", "leaprc.ff14SB", "", 0 );
unlink( "nab_tmp.pdb" );

return m;

};