// Program 8 - create a pseudoknot using distance bounds and refinement
molecule m;
float 	xyz[ dynamic ],f[ dynamic ],v[ dynamic ];
bounds 	b;
int	i, seqlen, ier;
float	fret;
string	seq, opt;

seq = "gcggaaacgccgcguaagcg";

seqlen = length(seq);

m = link_na("1", seq, "rna.amber94.rlb", "rna", "35");

allocate xyz[ 4*m.natoms ];
allocate f[ 4*m.natoms ];
allocate v[ 4*m.natoms ];

b = newbounds(m, "");

for ( i = 1; i <= seqlen; i = i + 1) {
	useboundsfrom(b, m, sprintf("1:%d:??,H?[^'T]", i), m,
		sprintf("1:%d:??,H?[^'T]", i), 0.0 );
}

setboundsfromdb(b, m, "1:1:", "1:2:", "arna.stack.db", 1.0);
setboundsfromdb(b, m, "1:2:", "1:3:", "arna.stack.db", 1.0);
setboundsfromdb(b, m, "1:3:", "1:18:", "arna.stack.db", 1.0);
setboundsfromdb(b, m, "1:18:", "1:19:", "arna.stack.db", 1.0);
setboundsfromdb(b, m, "1:19:", "1:20:", "arna.stack.db", 1.0);

setboundsfromdb(b, m, "1:8:", "1:9:", "arna.stack.db", 1.0);
setboundsfromdb(b, m, "1:9:", "1:10:", "arna.stack.db", 1.0);
setboundsfromdb(b, m, "1:10:", "1:11:", "arna.stack.db", 1.0);
setboundsfromdb(b, m, "1:11:", "1:12:", "arna.stack.db", 1.0);
setboundsfromdb(b, m, "1:12:", "1:13:", "arna.stack.db", 1.0);

setboundsfromdb(b, m, "1:1:", "1:13:", "arna.basepair.db", 1.0);
setboundsfromdb(b, m, "1:2:", "1:12:", "arna.basepair.db", 1.0);
setboundsfromdb(b, m, "1:3:", "1:11:", "arna.basepair.db", 1.0);

setboundsfromdb(b, m, "1:8:", "1:20:", "arna.basepair.db", 1.0);
setboundsfromdb(b, m, "1:9:", "1:19:", "arna.basepair.db", 1.0);
setboundsfromdb(b, m, "1:10:", "1:18:", "arna.basepair.db", 1.0);

tsmooth( b, 0.0005 );

opt = "seed=571, gdist=0, ntpr=50, k4d=2.0, randpair=5";
dg_options( b, opt );
embed(b, xyz );

ier = conjgrad( xyz, 4*m.natoms, fret, db_viol, 0.1, 10., 500 );

for ( i = 3000; i > 2800; i = i - 100 ){
	dg_options( b, "ntpr=1000, k4d=0.2" );
	mm_options( "ntpr_md=50, zerov=1, temp0=" + sprintf("%d.", i) );
	md( 4*m.natoms, 1000, xyz, f, v, db_viol );

	dg_options( b, "ntpr=1000, k4d=4.0" );
	mm_options( "zerov=0, temp0=0., tautp=0.3" );
	md( 4*m.natoms, 8000, xyz, f, v, db_viol );
}

setmol_from_xyzw( m, NULL, xyz );
putpdb( "rna_pseudoknot.pdb", m );