Skip to content

Commit

Permalink
fix ligand move distributions to be clearer
Browse files Browse the repository at this point in the history
  • Loading branch information
Justin Spiriti committed Jul 18, 2016
1 parent f704202 commit d5970af
Show file tree
Hide file tree
Showing 3 changed files with 34 additions and 32 deletions.
34 changes: 17 additions & 17 deletions src/init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -190,8 +190,8 @@ void simulation::process_commands(char * infname)
for (i=1; i<=NUM_MOVES; i++) {
prob[i]=0.0;
movesize[i]=0.0;
splitfrac[i]=0.0;
movesize2[i]=0.0;
large_dist_frac[i]=0.0;
movesize_large[i]=0.0;
}
for(;;) {
fgets(command,sizeof(command),input);
Expand All @@ -207,13 +207,13 @@ void simulation::process_commands(char * infname)
sscanf(command,"%s %lg %lg %lg %lg\n",word,&p,&size,&frac,&size2);
prob[move]=p;
movesize[move]=size;
splitfrac[move]=frac;
movesize2[move]=size2;
if (movesize2[move]<movesize[move]) {
temp=movesize2[move];
movesize2[move]=movesize[move];
large_dist_frac[move]=frac;
movesize_large[move]=size2;
if (movesize_large[move]<movesize[move]) {
temp=movesize_large[move];
movesize_large[move]=movesize[move];
movesize[move]=temp;
splitfrac[move]=1-splitfrac[move];
large_dist_frac[move]=1-large_dist_frac[move];
}
} else {
sscanf(command,"%s %lg %lg\n",word,&p,&size);
Expand All @@ -233,22 +233,22 @@ void simulation::process_commands(char * infname)
if ((move!=MOVE_LIGAND_TRANS) && (move!=MOVE_HEAVY_TRANS)) {
//strncpy(unit,"degrees\0",sizeof(unit));
movesize[move]*=DEG_TO_RAD;
movesize2[move]*=DEG_TO_RAD;
movesize_large[move]*=DEG_TO_RAD;
if (move==MOVE_LIGAND_ROT) {
printf("%.20s moves: Maximum size %.2f degrees Overall fraction %.2f%% Smaller distribution size %.2f degrees Split fraction %.2f%%\n",
mc_move_names[move],movesize2[move]*RAD_TO_DEG,prob[move]*100.0,movesize[move]*RAD_TO_DEG,splitfrac[move]*100.0);
printf("%.20s moves: Overall fraction %.2f%% -- Smaller max size %.2f degrees -- Larger max size %.2f degrees -- Large dist fraction %.2f%%\n",
mc_move_names[move],prob[move]*100.0,movesize[move]*RAD_TO_DEG,movesize_large[move]*RAD_TO_DEG,large_dist_frac[move]*100.0);
} else {
printf("%.20s moves: Maximum size %.2f degrees Overall fraction %.2f%%\n",
mc_move_names[move],movesize[move]*RAD_TO_DEG,prob[move]*100.0);
printf("%.20s moves: Overall fraction %.2f%% -- Maximum size %.2f degrees\n",
mc_move_names[move],prob[move]*100.0,movesize[move]*RAD_TO_DEG);
}
} else {
//strncpy(unit,"A\0",sizeof(unit);
if (move==MOVE_LIGAND_TRANS) {
printf("%.20s moves: Maximum size %.2f A Overall fraction %.2f%% Smaller distribution size %.2f A Split fraction %.2f%%\n",
mc_move_names[move],movesize2[move],prob[move]*100.0,movesize[move],splitfrac[move]*100.0);
printf("%.20s moves: Overall fraction %.2f%% -- Smaller max size %.2f A -- Larger max size %.2f A -- Large dist fraction %.2f%%\n",
mc_move_names[move],prob[move]*100.0,movesize[move],movesize_large[move],large_dist_frac[move]*100.0);
} else {
printf("%.20s moves: Maximum size %.2f degrees Overall fraction %.2f%%\n",
mc_move_names[move],movesize[move],prob[move]*100.0);
printf("%.20s moves: Overall fraction %.2f%% -- Maximum size %.2f A\n",
mc_move_names[move],prob[move]*100.0,movesize[move]);
}
}
}
Expand Down
8 changes: 4 additions & 4 deletions src/mc.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,8 +64,8 @@ class simulation {

long int nmcstep, nsave, nprint,ncheck,nprevstep; /*number of MC steps, frequency of saving conformations, frequency of printing*/
double prob[NUM_MOVES+1],cumprob[NUM_MOVES+1],movesize[NUM_MOVES+1];
//This allocates room for parameters for a split distribution. splitfrac is the fraction of t
double splitfrac[NUM_MOVES+1],movesize2[NUM_MOVES+1];
//This allocates room for parameters for a split distribution. large_dist_frac is the fraction of t
double large_dist_frac[NUM_MOVES+1],movesize_large[NUM_MOVES+1];
vector<mc_move> backbone_moves;
vector<mc_move> sidechain_moves;
vector<mc_move> backrub_moves;
Expand Down Expand Up @@ -116,8 +116,8 @@ class simulation {
void mcmove(int * movetype, subset * movedatoms, double * coords);
void rotate_atoms_by_axis(mc_move * move, const double angle, double * coords);
void rotate_atoms_by_point(subset atoms, const double * quat, const double * point, double * coords);
void do_ligand_trans(double movesize, double splitfrac, double movesize2, double * coords);
void do_ligand_rot(double movesize, double splitfrac, double movesize2, double * coords);
void do_ligand_trans(double movesize_small, double large_dist_frac, double movesize_large, double * coords);
void do_ligand_rot(double movesize_small, double large_dist_frac, double movesize_large, double * coords);
void heavy_atom_trans(subset * movedatoms, double movesize, double * coords);
void heavy_atom_rot(subset * movedatoms, double movesize, double * coords);
//dock prep related stuff -- see init.cpp
Expand Down
24 changes: 13 additions & 11 deletions src/mcmoves.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,25 +209,25 @@ void topology::generate_sidechain_moves(vector<mc_move> * sidechain_moves)
}
}

void simulation::do_ligand_trans(double movesize, double splitfrac, double movesize2, double * coords)
void simulation::do_ligand_trans(double movesize_small, double large_dist_frac, double movesize_large, double * coords)
{
double disp[3],m,r;
int k,iatom;
//construct a random vector within the unit sphere, and multiply by movesize
//to give a random displacement whose magnitude is no more than movesize.
r=genrand_real3();
if (r<splitfrac) {
rand_trans_vector(movesize,&disp[0]);
if (r<large_dist_frac) {
rand_trans_vector(movesize_large,&disp[0]);
} else {
rand_trans_vector(movesize2,&disp[0]);
rand_trans_vector(movesize_small,&disp[0]);
}
//translate all ligand atoms
for (iatom=0; iatom<top->natom; iatom++) if (top->ligand[iatom]) {
for (k=0; k<3; k++) coords[3*iatom+k]+=disp[k];
}
}

void simulation::do_ligand_rot(double movesize, double splitfrac, double movesize2, double * coords)
void simulation::do_ligand_rot(double movesize_small, double large_dist_frac, double movesize_large, double * coords)
{
double quat[4],com[3],mass,totmass,r;
int k, iatom;
Expand All @@ -241,13 +241,14 @@ void simulation::do_ligand_rot(double movesize, double splitfrac, double movesiz
}
for (k=0; k<3; k++) com[k]/=totmass;
r=genrand_real3();
if (r<splitfrac) {
rand_small_quat(movesize,&quat[0]);
if (r<large_dist_frac) {
rand_small_quat(movesize_large,&quat[0]);
} else {
rand_small_quat(movesize2,&quat[0]);
rand_small_quat(movesize_small,&quat[0]);
}
rotate_atoms_by_point(top->ligand,&quat[0],&com[0],coords);
}

//randomly move a heavy atom, togehter with any hydrogen atoms that may be bonded to it
void simulation::heavy_atom_trans(subset * movedatoms, double movesize, double * coords)
{
Expand All @@ -268,7 +269,8 @@ void simulation::heavy_atom_trans(subset * movedatoms, double movesize, double *
for (k=0; k<3; k++) coords[3*iatom+k]+=disp[k];
}
}
void simulation::heavy_atom_rot(subset * movedatoms, double movesize, double * coords)

void simulation::heavy_atom_rot(subset * movedatoms, double movesize, double * coords)
{
int iatom,iheavy,j,jatom,k;
double q[4];
Expand Down Expand Up @@ -314,11 +316,11 @@ void simulation::mcmove(int * movetype, subset * movedatoms, double * coords)
movecount=backrub_moves.size();
break;
case MOVE_LIGAND_TRANS:
do_ligand_trans(movesize[MOVE_LIGAND_TRANS],splitfrac[MOVE_LIGAND_TRANS],movesize2[MOVE_LIGAND_TRANS],coords);
do_ligand_trans(movesize[MOVE_LIGAND_TRANS],large_dist_frac[MOVE_LIGAND_TRANS],movesize_large[MOVE_LIGAND_TRANS],coords);
*movedatoms=top->ligand;
break;
case MOVE_LIGAND_ROT:
do_ligand_rot(movesize[MOVE_LIGAND_ROT],splitfrac[MOVE_LIGAND_ROT],movesize2[MOVE_LIGAND_ROT],coords);
do_ligand_rot(movesize[MOVE_LIGAND_ROT],large_dist_frac[MOVE_LIGAND_ROT],movesize_large[MOVE_LIGAND_ROT],coords);
*movedatoms=top->ligand;
break;
case MOVE_HEAVY_TRANS:
Expand Down

0 comments on commit d5970af

Please sign in to comment.