API tools faq
paste
Advertisement
Guest User

Nicolas

a guest
May 2nd, 2009
887
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 176.09 KB | None | 0 0
raw download clone embed print report
  1. /*
  2.  
  3.  $Id: main.cc,v 1.69 2007/05/01 02:40:12 garrett Exp $
  4.  
  5.  AutoDock
  6.  
  7.  Copyright (C) 1989-2007,  Garrett M. Morris, David S. Goodsell, Ruth Huey, Arthur J. Olson, Scott Halliday, Max Chang, Bill Hart, Richard Belew
  8.  All Rights Reserved.
  9.  
  10.  AutoDock is a Trade Mark of The Scripps Research Institute.
  11.  
  12.  This program is free software; you can redistribute it and/or
  13.  modify it under the terms of the GNU General Public License
  14.  as published by the Free Software Foundation; either version 2
  15.  of the License, or (at your option) any later version.
  16.  
  17.  This program is distributed in the hope that it will be useful,
  18.  but WITHOUT ANY WARRANTY; without even the implied warranty of
  19.  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20.  GNU General Public License for more details.
  21.  
  22.  You should have received a copy of the GNU General Public License
  23.  along with this program; if not, write to the Free Software
  24.  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  25.  
  26.  */
  27.  
  28. #ifdef HAVE_CONFIG_H
  29. #include <config.h>
  30. #endif
  31.  
  32. #include <math.h>
  33. #include <sys/types.h> // time_t time(time_t *tloc);
  34. #include <time.h>      // time_t time(time_t *tloc);
  35. #include <sys/times.h>
  36. #include <stdlib.h>
  37. #include <string.h>
  38. #include <sys/param.h>
  39. #include <ctype.h> // tolower
  40. #include <unistd.h> // sysconf
  41.  
  42. /* the BOINC API header file */
  43. #ifdef BOINC
  44. #include "diagnostics.h"
  45. #include "boinc_api.h"
  46. #include "filesys.h"                 // boinc_fopen(), etc... */
  47. #endif
  48.  
  49. #include "coliny.h"
  50. #include "hybrids.h"
  51. #include "ranlib.h"
  52. #include "gs.h"
  53. #include "ls.h"
  54. #include "rep.h"
  55. #include "support.h"
  56. #include "distdepdiel.h"
  57. #include "calculateEnergies.h"
  58. #include "conformation_sampler.h"
  59. #include "main.h"
  60.  
  61. extern int debug;
  62. extern int keepresnum;
  63. extern Real idct;
  64. extern Eval evaluate;
  65. extern Linear_FE_Model AD4;
  66. extern Real nb_group_energy[3]; ///< total energy of each nonbond group (intra-ligand, inter, and intra-receptor)
  67. extern int Nnb_array[3];  ///< number of nonbonds in the ligand, intermolecular and receptor groups
  68.  
  69. static const char* const ident[] = {ident[1], "@(#)$Id: main.cc,v 1.69 2007/05/01 02:40:12 garrett Exp $"};
  70.  
  71.  
  72. int sel_prop_count = 0;
  73.  
  74.  
  75. int main (int argc, char * const argv[], char * const envp[])
  76.  
  77. /*******************************************************************************
  78. **      Name: main  (AutoDock)                                                **
  79. **  Function: Performs Automated Docking of Small Molecule into Macromolecule **
  80. ** Copyright: (C) 1994-2005 TSRI, Arthur J. Olson's Labortatory.              **
  81. **____________________________________________________________________________**
  82. **   Authors: Garrett Matthew Morris, Current C/C++ version 4.0               **
  83. **                                       e-mail: garrett@scripps.edu          **
  84. **                                                                            **
  85. **            David Goodsell, Orignal FORTRAN version 1.0                     **
  86. **                                       e-mail: goodsell@scripps.edu         **
  87. **                                                                            **
  88. **            The Scripps Research Institute                                  **
  89. **            Department of Molecular Biology, MB5                            **
  90. **            10550 North Torrey Pines Road                                   **
  91. **            La Jolla, CA 92037.                                             **
  92. **                                                                            **
  93. **      Date: 02/10/04                                                        **
  94. **____________________________________________________________________________**
  95. **    Inputs: Control file, Small Molecule PDBQT file, macromolecular grid map**
  96. **            files.                                                          **
  97. **   Returns: Autodock Log File, includes docked conformation clusters (PDBQT)**
  98. **   Globals: X, Y, Z, SPACE, MAX_GRID_PTS, NEINT, MAX_ATOMS,                 **
  99. **            MAX_MAPS, ATOM_MAPS, A_DIVISOR, LINE_LEN,TRUE,FALSE             **
  100. **            programname, command_mode,                                      **
  101. **            command_in_fp, command_out_fp, parFile, logFile.                **
  102. **____________________________________________________________________________**
  103. ** Modification Record                                                        **
  104. ** Date     Inits   Comments                                                  **
  105. ** 09/06/95 RSH     Added code to handle GA/LS stuff                          **
  106. **          DSG     Quaternion rotations                                      **
  107. **          DSG     Generates torsions from annotated pdb list                **
  108. **          DSG     Generates internal energies                               **
  109. **          DSG     Performs a limited Cluster analysis of conformations      **
  110. ** 05/07/92 GMM     C translation                                             **
  111. ** 05/14/92 GMM     Time-dependent seed in random-number generation           **
  112. ** 10/29/92 GMM     Application Visualization System (AVS) readable grid      **
  113. **                  display file input.                                       **
  114. **                  [AVS is a trademark of Stardent Computer Inc.]            **
  115. **                  Also added the 'total_charge' check.                      **
  116. ** 11/19/93 GMM     #ifdef NOSQRT, with non-square-rooting acceleration.      **
  117. ** 09/26/94 GMM     Cluster analysis now outputs RMS deviations.              **
  118. ** 09/28/94 GMM     Modularized code.                                         **
  119. ** 10/02/94 GMM     Distance constraints added, for Ed Moret. Accelerated.    **
  120. ** 09/06/95 RSH     Incorporation of GA/SW tokens                             **
  121. *******************************************************************************/
  122.  
  123. {
  124. //   ATOM_MAPS
  125. //
  126. char            atm_typ_str[ATOM_MAPS]; //  "atm_typ_str" (in AD3) used to serve a similar role to "atom_type_name" (in AD4).
  127.  
  128. //   MAX_MAPS
  129. //
  130. char            *ligand_atom_type_ptrs[MAX_MAPS]; /* array of ptrs used to parse input line of atom type names */
  131. ParameterEntry  parameterArray[MAX_MAPS];
  132.  
  133. //   MAX_GRID_PTS & MAX_MAPS
  134. //
  135. static Real map[MAX_GRID_PTS][MAX_GRID_PTS][MAX_GRID_PTS][MAX_MAPS];
  136. // double *map;  // Use this with malloc...
  137.  
  138. GridMapSetInfo *info;  // this information is from the AVS field file
  139.  
  140. //   MAX_ATOMS
  141. //
  142. char atomstuff[MAX_ATOMS][MAX_CHARS];
  143. char pdbaname[MAX_ATOMS][5];
  144. Real crdpdb[MAX_ATOMS][SPACE];  // original PDB coordinates from input
  145. Real crdreo[MAX_ATOMS][SPACE];  // reoriented coordinates
  146. Real crd[MAX_ATOMS][SPACE];     // current coordinates
  147. Real charge[MAX_ATOMS];
  148. Real abs_charge[MAX_ATOMS];
  149. Real qsp_abs_charge[MAX_ATOMS];
  150. Real elec[MAX_ATOMS];
  151. Real emap[MAX_ATOMS];
  152. int type[MAX_ATOMS];
  153. int bond_index[MAX_ATOMS];
  154. int ignore_inter[MAX_ATOMS];
  155. Atom atoms[MAX_ATOMS];
  156.  
  157. //   MAX_TORS
  158. //
  159. int  tlist[MAX_TORS][MAX_ATOMS];
  160. Real vt[MAX_TORS][SPACE];
  161. Real F_TorConRange[MAX_TORS][MAX_TOR_CON][2];
  162. unsigned short  US_TorE[MAX_TORS];
  163. Boole B_isTorConstrained[MAX_TORS];
  164. int N_con[MAX_TORS];
  165. unsigned short US_torProfile[MAX_TORS][NTORDIVS];
  166.  
  167. //   MAX_NONBONDS
  168. //
  169. NonbondParam *nonbondlist = new NonbondParam[MAX_NONBONDS];
  170.  
  171. //   LINE_LEN
  172. //
  173. char error_message[LINE_LEN];
  174. char message[LINE_LEN];
  175. char line[LINE_LEN];
  176. char torfmt[LINE_LEN];
  177.  
  178. //   MAX_CHARS
  179. //
  180. char FN_clus[MAX_CHARS];
  181. char FN_watch[MAX_CHARS];
  182. //char FN_gdfld[MAX_CHARS]  // now part of the GridMapSetInfo structure;
  183. //char FN_gpf[MAX_CHARS];  // now part of the GridMapSetInfo structure
  184. char FN_ligand[MAX_CHARS];
  185. char dummy_FN_ligand[MAX_CHARS];
  186. char FN_trj[MAX_CHARS];
  187. //char FN_receptor[MAX_CHARS];// now part of the GridMapSetInfo structure
  188. char FN_rms_ref_crds[MAX_CHARS];
  189. char FN_parameter_library[MAX_CHARS];
  190. char hostnm[MAX_CHARS];
  191. char param[2][MAX_CHARS];
  192. char c_mode_str[MAX_CHARS];
  193. char FN_pop_file[MAX_CHARS];
  194. char FN_flexres[MAX_CHARS];
  195. char rms_atoms_cmd[MAX_CHARS];
  196. char confsampler_type[MAX_CHARS];
  197.  
  198. //   MAX_RECORDS
  199. //
  200. char PDBQT_record[MAX_RECORDS][LINE_LEN];
  201.  
  202. //   SPACE
  203. //
  204. Real lig_center[SPACE];
  205. //Real map_center[SPACE];
  206.  
  207. //   MAX_RUNS
  208. //
  209. Real econf[MAX_RUNS];  // this is the list of energies printed in the histogram in "analysis"
  210. State sHist[MAX_RUNS];  /* qtnHist[MAX_RUNS][QUAT],torHist[MAX_RUNS][MAX_TORS];*/
  211. State sUnbound; // State of the unbound ligand's conformation
  212. State sUnbound_ext; // State of the unbound ligand's conformation after extended-conformation search
  213. //// UNCOMMENT if using Step 2 in unbound calculation ---> State sUnbound_ls; // State of the unbound ligand's conformation after a local search
  214. State sUnbound_ad; // State of the unbound ligand's conformation after an AutoDock search
  215.  
  216. char out_acc_rej = '?';
  217. char timeSeedIsSet[2];
  218. char selminpar = 'm';
  219. char S_contype[8];
  220.  
  221. static ParameterEntry * foundParameter;
  222.  
  223. Real cA;
  224. Real cB;
  225. Real cA_unbound = 392586.8;  // repulsive
  226. Real cB_unbound = 0.0; // attractive
  227. Real epsij;
  228. Real F_A;
  229. Real F_Aova;
  230. Real F_tor;
  231. Real F_torPref;
  232. Real F_torHWdth;
  233. Real Rij;
  234. Real sqlower;
  235. Real squpper;
  236. Real tmpconst;
  237.  
  238. // Distance-dependence in Desolvation Term
  239. const double sigma = 3.6L;
  240. const double qsolpar = 0.01097L;
  241.  
  242. // ELECSCALE converts between CGS units and SI units;
  243. // see, e.g. p 254, "Molecular Modeling and Simulation", by Tamar Schlick, Springer.
  244. //
  245. // Units of ELECSCALE are (Kcal/mol ) * (Angstrom / esu^2)
  246. // and this allows us to use distances in  Angstroms and charges in esu...
  247. const Real ELECSCALE = 332.06363;  
  248.  
  249. // const Real ELECSCALE = 83.0159075;   this ELECSCALE (corresponding to eps(r) = 1/4r) gives -7.13 kcal/mol for 1pgp Tests/test_autodock4.py
  250.  
  251. // i
  252. double Ri, epsi, Ri_hb, epsi_hb;
  253. hbond_type hbondi;
  254.  
  255. // j
  256. double Rj, epsj, Rj_hb, epsj_hb;
  257. hbond_type hbondj;
  258.  
  259. Real scale_1_4 = 0.5;
  260. Real c=0.0;
  261. Real clus_rms_tol = 0.0;
  262. Real e0max = BIG;
  263. Real eintra = 0.0;  // sum of intramolecular energy for the ligand plus that of the protein
  264. Real einter = 0.0; // intermolecular energy between the ligand and the protein
  265. Real etotal = 0.0;
  266. Real AD3_FE_coeff_estat   = 1.000;
  267. Real qtwFac = 1.0;
  268. Real qtwStep0 = DegreesToRadians( 5.0 );
  269. Real qtwStepFinal = DegreesToRadians( 0.5 );
  270. Real maxrad = -1.0;
  271. Real r2sum=0.0;
  272. // Real RJ = 8.31441;     // in J/K/mol, Gas Constant, Atkins Phys.Chem., 2/e
  273. // Real Rcal = 1.9871917; // in cal/K/mol, Gas Constant, RJ/4.184
  274. // Real T0K = 273.15;        // 0 degrees Celsius, in K
  275. Real RTreduc = 1.0;
  276. // Real spacing = 0.0;// now part of the GridMapSetInfo structure
  277. Real RT0 = 616.0;
  278. Real RTFac = 0.95;
  279. Real torsdoffac = 0.3113;
  280. Real torsFreeEnergy = 0.0;
  281. Real torFac = 1.0;
  282. Real torStep0 = 5.0;
  283. Real torStepFinal = 5.0;
  284. Real trnFac = 1.0;
  285. Real trnStep0 = 0.2;
  286. Real trnStepFinal = 0.2;
  287. Real WallEnergy = 1.0e8; /* Energy barrier beyond walls of gridmaps. */
  288. //  The GA Stuff
  289. Real m_rate = 0.02;
  290. Real c_rate = 0.80;
  291. Real alpha = 0;
  292. Real beta = 1;
  293. Real search_freq = 0.06;
  294. Real rho = 1.0;
  295. Real lb_rho = 0.01;
  296. Real *rho_ptr = NULL;
  297. Real *lb_rho_ptr = NULL;
  298. Real unbound_internal_FE = 0.0;
  299. Real unbound_ext_internal_FE = 0.0;
  300. Real unbound_ad_internal_FE = 0.0;
  301. Real unbound_internal_FE_saved = 0.0;
  302. Real emap_total = 0.;
  303. Real elec_total = 0.;
  304. Real charge_total = 0.;
  305. Real etot = 0.;
  306.  
  307. EnergyBreakdown eb;
  308.  
  309. initialise_energy_breakdown(&eb, 0, 0);
  310.  
  311. unsigned int outputEveryNgens = 100;
  312.  
  313. Boole B_atom_types_found = FALSE;
  314. Boole B_isGaussTorCon = FALSE;
  315. Boole B_constrain_dist;
  316. Boole B_either = FALSE;
  317. Boole B_calcIntElec = FALSE;
  318. Boole B_calcIntElec_saved = FALSE;
  319. Boole B_write_trj = FALSE;
  320. Boole B_watch = FALSE;
  321. Boole B_acconly = FALSE;
  322. Boole B_cluster_mode = FALSE;
  323. Boole B_havemap = FALSE;
  324. Boole B_havenbp = FALSE;
  325. Boole B_haveCharges;
  326. Boole B_linear_schedule = FALSE;
  327. Boole B_qtwReduc = FALSE;
  328. Boole B_selectmin = FALSE;
  329. Boole B_symmetry_flag = TRUE;
  330. Boole B_tempChange = TRUE;
  331. Boole B_torReduc = FALSE;
  332. Boole B_trnReduc = FALSE;
  333. Boole B_write_all_clusmem = FALSE;
  334. Boole B_ShowTorE = FALSE;
  335. Boole B_RandomTran0 = FALSE;
  336. Boole B_RandomQuat0 = FALSE;
  337. Boole B_RandomDihe0 = FALSE;
  338. Boole B_CalcTrnRF = FALSE;
  339. Boole B_CalcQtwRF = FALSE;
  340. Boole B_CalcTorRF = FALSE;
  341. Boole B_charMap = FALSE;
  342. Boole B_include_1_4_interactions = FALSE;  // This was the default behaviour in previous AutoDock versions (1 to 3).
  343. Boole B_found_move_keyword = FALSE;
  344. Boole B_found_ligand_types = FALSE;
  345. Boole B_found_elecmap = FALSE;
  346. Boole B_found_desolvmap = FALSE;
  347. Boole B_use_non_bond_cutoff = TRUE;
  348. Boole B_have_flexible_residues = FALSE;  // if the receptor has flexible residues, this will be set to TRUE
  349. Boole B_rms_atoms_ligand_only = TRUE;  // cluster on the ligand atoms only
  350. Boole B_reorient_random = FALSE; // if true, create a new random orientation before docking
  351. int atm1=0;
  352. int atm2=0;
  353. int a1=0;
  354. int a2=0;
  355. int atomC1;
  356. int atomC2;
  357. int dpf_keyword = -1;
  358. //int gridpts1[SPACE];  // now part of the GridMapSetInfo structure
  359. //int gridpts[SPACE];  // now part of the GridMapSetInfo structure
  360. int Htype = 0;
  361. int ncycles = -1;
  362. int iCon=0;
  363. int indcom = 0;
  364. int ligand_is_inhibitor = 1;
  365. int ltorfmt = 4;
  366. int nruns = 0;
  367. int nstepmax = -1;
  368. int naccmax = 0;
  369. int natom = 0;
  370. int nconf = 0;
  371. int ncycm1 = 1;
  372. int ndihed = 0;
  373. int nlig = 0;
  374. int nres = 0;
  375. int nmol = 0;
  376. int Nnb = 0;
  377. int nrejmax = 0;
  378. int ntor = 0;
  379. int ntor1 = 1;
  380. int ntor_ligand = 0;
  381. int ntorsdof = 0;
  382. int num_maps = 0;
  383. int num_atom_types = 0;
  384. int nval = 0;
  385. int outlev = -1;
  386. int retval = 0;
  387. int status = 0;
  388. int trj_end_cyc = 0;
  389. int trj_begin_cyc = 0;
  390. int trj_freq = 0;
  391. int xA = 12;
  392. int xB = 6;
  393. int xA_unbound = 12;
  394. int xB_unbound = 6;
  395. int I_tor;
  396. int I_torBarrier;
  397. int MaxRetries = 1000; /* Default maximum number of retries for ligand init. */
  398. int OutputEveryNTests = 1000;
  399. int NumLocalTests = 10;
  400. int maxTests = 10000;
  401. int parameter_library_found = 0;
  402. /* int beg; */
  403. /* int end; */
  404. /* int imol = 0; */
  405. int outside = FALSE;
  406. int atoms_outside = FALSE;
  407. // unsigned int min_evals_unbound =  250000;
  408. unsigned int max_evals_unbound = 1000000;
  409. int saved_sInit_ntor = 0;
  410. int confsampler_samples = 0;
  411.  
  412. unsigned short US_energy;
  413. unsigned short US_tD;
  414. unsigned short US_torBarrier = TORBARMAX;
  415. unsigned short US_min = TORBARMAX;
  416.  
  417. register int i = 0;
  418. register int j = 0;
  419. register int k = 0;
  420. //register int m = 0;
  421. register int xyz = 0;
  422. int j1 = 1;
  423.  
  424. State sInit;            /* Real qtn0[QUAT], tor0[MAX_TORS]; */
  425.  
  426. Quat q_reorient;
  427.  
  428. Molecule ligand;        /* ligand */
  429.  
  430. static Real F_A_from;
  431. static Real F_A_to;
  432. static Real F_lnH;
  433. static Real F_W;
  434. static Real F_hW;
  435. static FourByteLong clktck = 0;
  436.  
  437. static Real version_num = 4.00;
  438.  
  439. struct tms tms_jobStart;
  440. struct tms tms_gaStart;
  441. struct tms tms_gaEnd;
  442.  
  443. Clock  jobStart;
  444. Clock  gaStart;
  445. Clock  gaEnd;
  446.  
  447. time_t time_seed;
  448.  
  449. EnergyTables *ad_energy_tables;  // Holds vdw+Hb, desolvation & dielectric lookup tables
  450. EnergyTables *unbound_energy_tables;  // Use for computing unbound energy & conformation
  451.  
  452. Statistics map_stats;
  453.  
  454. //  The GA Stuff
  455. FourByteLong seed[2];
  456. unsigned int pop_size = 200;
  457. unsigned int num_generations = 0;  //  Don't terminate on the basis of number of generations
  458. unsigned int num_evals = 250000;
  459. unsigned int num_evals_unbound = num_evals;
  460. unsigned int max_its = 300;
  461. unsigned int max_succ = 4;
  462. unsigned int max_fail = 4;
  463. int window_size = 10;
  464. int low = 0;
  465. int high = 100;
  466. int elitism = 1;
  467.  
  468.  
  469. Selection_Mode s_mode = Proportional;
  470. Xover_Mode c_mode = TwoPt;  //  can be: OnePt, TwoPt, Uniform or Arithmetic
  471. Worst_Mode w_mode = AverageOfN;
  472. EvalMode e_mode = Normal_Eval;
  473. Global_Search *GlobalSearchMethod = NULL;
  474. Local_Search *LocalSearchMethod = NULL;
  475. //// Local_Search *UnboundLocalSearchMethod = NULL;
  476.  
  477. info = (GridMapSetInfo *) malloc( sizeof(GridMapSetInfo) );
  478. ad_energy_tables = (EnergyTables *) malloc( sizeof(EnergyTables) );
  479. unbound_energy_tables = (EnergyTables *) malloc( sizeof(EnergyTables) );
  480.  
  481. // Create a coordinate at the origin:
  482. Coord origin;
  483. origin.x = 0.;
  484. origin.y = 0.;
  485. origin.z = 0.;
  486.  
  487. //______________________________________________________________________________
  488. /*
  489. ** Get the time at the start of the run...
  490. */
  491.  
  492. jobStart = times( &tms_jobStart );
  493.  
  494.  
  495. //_____________________________________________________________________________
  496. /*
  497. ** Boinc initialization
  498. */
  499. #ifdef BOINC
  500.     int flags = 0;
  501.     int rc;
  502.     flags =
  503.       BOINC_DIAG_DUMPCALLSTACKENABLED |
  504.       BOINC_DIAG_HEAPCHECKENABLED |
  505.       BOINC_DIAG_REDIRECTSTDERR |
  506.       BOINC_DIAG_REDIRECTSTDOUT ;
  507.     boinc_init_diagnostics(flags);
  508.  
  509. #ifdef BOINCCOMPOUND
  510.     BOINC_OPTIONS options;
  511.     options.main_program = false;
  512.     options.check_heartbeat = false; // monitor does check heartbeat
  513.     options.handle_trickle_ups = false;
  514.     options.handle_trickle_downs = false;
  515.     options.handle_process_control = false;
  516.     options.send_status_msgs = true;// only the worker programs (i.e. model) sends status msgs
  517.     options.direct_process_action = true;// monitor handles suspend/quit, but app/model doesn't
  518.     // Initialisation of Boinc
  519.     rc =  boinc_init_options(options); //return 0 for success
  520.     if( rc ){
  521.       fprintf(stderr,"BOINC_ERROR: boinc_init_options() failed \n");
  522.       exit(rc);
  523.     }
  524.  
  525. #else
  526.     // All BOINC applications must initialize the BOINC interface:
  527.     rc = boinc_init();
  528.     if (rc){
  529.       fprintf(stderr, "BOINC_ERROR: boinc_init() failed.\n");
  530.       exit(rc);
  531.     }
  532. #endif
  533. #endif
  534.  
  535.  
  536. //______________________________________________________________________________
  537. /*
  538. ** Parse the arguments in the command line...
  539. */
  540.  
  541. if ( setflags(argc,argv) == -1) {
  542.     exit(-1);
  543. } /* END PROGRAM */
  544.  
  545.  
  546. //______________________________________________________________________________
  547. /*
  548. ** Initialize torsion arrays and constants.
  549. */
  550.  
  551. (void) strcpy( torfmt, "%*s\0" ); /* len(torfmt) is 4 chars */
  552.  
  553. for (j = 0;  j < MAX_ATOMS;  j++ ) {
  554.     type[j] = 0;
  555.     ignore_inter[j] = 0;
  556. }
  557.  
  558. for (i = 0; i  < MAX_TORS;  i++ ) {
  559.     for (j = 0;  j < MAX_ATOMS;  j++ ) {
  560.         tlist[i][j] = 0;
  561.     }
  562. }
  563.  
  564. for (i = 0; i  < MAX_TORS;  i++ ) {
  565.     B_isTorConstrained[i] = 0;
  566.     US_torProfile[i][0] = 0;
  567.     N_con[i] = 0;
  568. }
  569.  
  570. initialiseState( &sInit );
  571. initialiseState( &(ligand.S) );
  572.  
  573. initialiseQuat( &q_reorient );
  574.  
  575. F_W      =  360.0 / NTORDIVS;
  576. F_hW     =  F_W  / 2.0;
  577. F_A_from = -360.0 + F_hW;
  578. F_A_to   =  360.0 + F_hW;
  579.  
  580. for (k = 0; k < MAX_RUNS; k++) {
  581.     for (i = 0; i  < MAX_TORS;  i++ ) {
  582.         sHist[k].tor[i] = 0.0;
  583.     }
  584. }
  585.  
  586. for (i = 0; i < MAX_TORS;  i++ ) {
  587.     if ( (ltorfmt += 4) > LINE_LEN ) {
  588.         prStr( error_message, "%s:  ERROR: MAX_TORS = %d torsions declared in \"constants.h\";\n\t LINE_LEN = %d, Therefore you must change \"LINE_LEN\" to exceed %d...\n", programname, MAX_TORS, LINE_LEN, 4+4*MAX_TORS );
  589.         stop( error_message );
  590.         exit( -1 );
  591.     } else {
  592.         (void) strcat( torfmt, " %lf\0" );  /* add on 4 chars  for each new torsion... */
  593.     }
  594. } /* len(torfmt) is 4+4*MAX_TORS chars */
  595.  
  596. for (j = 0; j < MAX_NONBONDS; j++) {
  597.     nonbondlist[j].a1 = nonbondlist[j].a2 = 0;
  598. }
  599.  
  600. for (j = 0; j < MAX_RUNS; j++) {
  601.     // isort[j] = j;
  602.     econf[j] = 0.0;
  603. }
  604.  
  605. B_constrain_dist = B_haveCharges = FALSE;
  606. ntor1 = ntor = atomC1 = atomC2 = 0;
  607. sqlower = squpper = 0.0;
  608.  
  609. timeSeedIsSet[0] = 'F';
  610. timeSeedIsSet[1] = 'F';
  611.  
  612. if (clktck == 0) {        /* fetch clock ticks per second first time */
  613.     if ( (clktck = sysconf(_SC_CLK_TCK)) < (FourByteLong)0L) {
  614.         stop("\"sysconf(_SC_CLK_TCK)\" command failed in \"main.c\"\n");
  615.         exit( -1 );
  616.     } else {
  617.         idct = (Real)1.0 / (Real)clktck;
  618.         if (debug) {
  619.             pr(logFile, "N.B. debug is on and set to %d\n\n", debug);
  620.             pr(logFile, "\n\nFYI:  Number of clock ticks per second = %d\n", (int)clktck);
  621.             pr(logFile, "FYI:  Elapsed time per clock tick = %.3e milli-seconds\n\n\n\n", idct * 1000. );
  622.         }
  623.     }
  624. }
  625.  
  626. (void) strcpy(FN_rms_ref_crds,"unspecified filename\0");
  627.  
  628.  
  629. //______________________________________________________________________________
  630. /*
  631. ** log(x): compute the natural (base e) logarithm of x,
  632. */
  633.  
  634. F_lnH = ((Real)log(0.5));
  635.  
  636. //______________________________________________________________________________
  637. /*
  638. ** Determine output level before we ouput anything.  
  639. ** We must parse the entire DPF -- silently -- for any outlev settings.
  640. */
  641.  
  642. while( fgets(line, LINE_LEN, parFile) != NULL ) { /* PARSING-DPF parFile */
  643.     dpf_keyword = parse_dpf_line( line );
  644.  
  645.     switch( dpf_keyword ) {
  646.     case DPF_OUTLEV:
  647.         /*
  648.         **  outlev
  649.         **  Output level,
  650.         */
  651.         retval = sscanf( line, "%*s %d", &outlev );
  652.         switch ( outlev ) {
  653.         case -1:
  654.             outputEveryNgens = (unsigned int) OUTLEV0_GENS;
  655.             break;
  656.         case 0:
  657.             outputEveryNgens = (unsigned int) OUTLEV0_GENS;
  658.             break;
  659.         case 1:
  660.             outputEveryNgens = (unsigned int) OUTLEV1_GENS;
  661.             break;
  662.         case 2:
  663.         default:
  664.             outputEveryNgens = (unsigned int) OUTLEV2_GENS;
  665.             break;
  666.         }
  667.         break;
  668.  
  669.     case DPF_FLEXRES:
  670.         // The DPF specifies a flexible residues file
  671.         // -- set a flag
  672.         // -- get the filename
  673.         B_have_flexible_residues = TRUE;
  674.         (void) sscanf( line, "%*s %s", FN_flexres );
  675.         break;
  676.  
  677.     default:
  678.         break;
  679.     } // switch( dpf_keyword )
  680. } // while
  681.  
  682. // Rewind DPF, so we can resume normal parsing
  683. (void) rewind( parFile );
  684.  
  685. //______________________________________________________________________________
  686. /*
  687. ** Output banner...
  688. */
  689.  
  690. banner( version_num );
  691.  
  692. (void) fprintf(logFile, "                           $Revision: 1.69 $\n\n");
  693. (void) fprintf(logFile, "                   Compiled on %s at %s\n\n\n", __DATE__, __TIME__);
  694.  
  695.  
  696. //______________________________________________________________________________
  697. /*
  698. ** Print the time and date when the file was created...
  699. */
  700.  
  701. pr( logFile, "This file was created at:\t\t\t" );
  702. printdate( logFile, 1 );
  703.  
  704. (void) strcpy(hostnm, "unknown host\0");
  705.  
  706. if (gethostname( hostnm, MAX_CHARS ) == 0) {
  707.     pr( logFile, "                   using:\t\t\t\"%s\"\n", hostnm );
  708. }
  709.  
  710. pr( logFile, "\nNOTE: \"rus\" stands for:\n\n      r = Real, wall-clock or elapsed time;\n      u = User or cpu-usage time;\n      s = System time\n\nAll timings are in seconds, unless otherwise stated.\n\n\n" );
  711.  
  712. //______________________________________________________________________________
  713.  
  714. (void) fprintf(logFile, "\n      ________________________________________________________________\n\n");
  715. (void) fprintf(logFile, "                   SETTING UP DEFAULT PARAMETER LIBRARY\n");
  716. (void) fprintf(logFile, "      ________________________________________________________________\n\n\n");
  717.  
  718. //______________________________________________________________________________
  719. //
  720. // Read in default parameters
  721. //
  722. setup_parameter_library(outlev);
  723.  
  724. //
  725. // Compute the look-up table for the distance-dependent dielectric function
  726. //
  727. (void) fprintf(logFile, "\n\nPreparing Energy Tables for Bound Calculation:\n\n");
  728. setup_distdepdiel(outlev, ad_energy_tables);
  729. (void) fprintf(logFile, "Preparing Energy Tables for Unbound Calculation:\n\n");
  730. setup_distdepdiel(outlev, unbound_energy_tables);
  731.  
  732. //______________________________________________________________________________
  733.  
  734. (void) fprintf(logFile, "\n      ___________________________________________________\n\n");
  735. (void) fprintf(logFile,   "             PARSING INPUT DOCKING PARAMETER FILE\n");
  736. (void) fprintf(logFile,   "      ___________________________________________________\n\n");
  737.  
  738. //______________________________________________________________________________
  739. /*
  740. ** (Note: "dock_param_fn" set in "setflags.c"...)
  741. */
  742. pr( logFile, "Docking parameter file (DPF) used for this docking:\t\t%s\n\n", dock_param_fn );
  743.  
  744. //______________________________________________________________________________
  745. /*
  746. ** Start reading in the DPF parameter/run-control file,
  747. */
  748.  
  749. while( fgets(line, LINE_LEN, parFile) != NULL ) { /* PARSING-DPF parFile */
  750.     // "line" is a string containing the current line of the input DPF.
  751.  
  752.     dpf_keyword = parse_dpf_line( line );
  753.  
  754.     switch( dpf_keyword ) {
  755.         case -1:
  756.             pr( logFile, "DPF> %s", line );
  757.             prStr( error_message,"%s: WARNING: Unrecognized keyword in docking parameter file.\n", programname );
  758.             pr_2x( logFile, stderr, error_message );
  759.             continue;
  760.             /* break; */
  761.  
  762.         case DPF_NULL:
  763.         case DPF_COMMENT:
  764.             pr( logFile, "DPF> %s", line );
  765.             (void) fflush(logFile);
  766.             break;
  767.  
  768.         default:
  769.             pr( logFile, "\n\nDPF> %s\n", line );
  770.             indcom = strindex( line, "#" );
  771.             if (indcom != -1) {
  772.                 line[ indcom ] = '\0'; /* Truncate "line" at the comment */
  773.             }
  774.             (void) fflush(logFile);
  775.             break;
  776.     } /* switch */
  777.  
  778.     switch( dpf_keyword ) {
  779.  
  780. //______________________________________________________________________________
  781.  
  782.     case DPF_NULL:
  783.     case DPF_COMMENT:
  784.         break;
  785.  
  786. //______________________________________________________________________________
  787.  
  788.     case DPF_OUTLEV:
  789.         /*
  790.         **  outlev
  791.         **  Output level,
  792.         */
  793.         retval = sscanf( line, "%*s %d", &outlev );
  794.         switch ( outlev ) {
  795.         case -1:
  796.             pr( logFile, "Output Level = -1.  ONLY STATE VARIABLES OUTPUT, NO COORDINATES.\n" );
  797.             outputEveryNgens = (unsigned int) OUTLEV0_GENS;
  798.             break;
  799.         case 0:
  800.             pr( logFile, "Output Level = 0.\n" );
  801.             outputEveryNgens = (unsigned int) OUTLEV0_GENS;
  802.             break;
  803.         case 1:
  804.             pr( logFile, "Output Level = 1.  MINIMUM OUTPUT DURING DOCKING.\n" );
  805.             outputEveryNgens = (unsigned int) OUTLEV1_GENS;
  806.             break;
  807.         case 2:
  808.         default:
  809.             pr( logFile, "Output Level = 2.  FULL OUTPUT DURING DOCKING.\n" );
  810.             outputEveryNgens = (unsigned int) OUTLEV2_GENS;
  811.             break;
  812.         }
  813.         pr( logFile, "\n\tOutput every %u generations.\n", outputEveryNgens );
  814.         (void) fflush(logFile);
  815.         break;
  816.  
  817. /*____________________________________________________________________________*/
  818.  
  819.     case DPF_PARAMETER_LIBRARY:
  820.         /*
  821.         ** parameter_file AD4_parameters.dat
  822.         **  or
  823.         ** parameter_library AD4_parameters.dat
  824.         **
  825.         ** initial implementation based on hsearch was suggested by Mike Pique
  826.         */
  827.  
  828.         parameter_library_found = sscanf( line, "%*s %s", FN_parameter_library );
  829.         (void) fflush(logFile);
  830.  
  831.         read_parameter_library(FN_parameter_library, outlev);
  832.  
  833.         break;
  834.  
  835. /*____________________________________________________________________________*/
  836.  
  837.     case DPF_INCLUDE_1_4_INTERACTIONS:
  838.         /*
  839.          * include_1_4_interactions 0.5
  840.          *
  841.          * Set the Boolean variable, B_include_1_4_interactions, to TRUE.
  842.          *
  843.          * NOTE:  You must use this command _before_ the "move ligand.pdbqt"
  844.          *        command, since "include_1_4_interactions" affects how the Ligand
  845.          *        PDBQT specified by the "move" command will be interpreted.
  846.          */
  847.         if (B_found_move_keyword == TRUE) {
  848.             // If we have found the move keyword already, warn the user
  849.             // that this command ("include_1_4_interactions 0.5") should have
  850.             // been given before this!
  851.             pr(logFile, "\nWARNING:  This command will be ignored.\n\nYou must put this command _before_ the \"move ligand.pdbqt\" command, since this command affects how the PDBQT file will be interpreted.\n\n");
  852.         }
  853.         (void) sscanf( line, "%*s " FDFMT, &scale_1_4 );
  854.         B_include_1_4_interactions = TRUE;
  855.         print_1_4_message(logFile, B_include_1_4_interactions, scale_1_4);
  856.         break;
  857.  
  858. //______________________________________________________________________________
  859.  
  860.     case DPF_INTELEC:
  861.         /*
  862.         **  intelec
  863.         **  Calculate internal electrostatic energies...
  864.         */
  865.         B_calcIntElec = TRUE;
  866.         if (outlev >= 0) {
  867.             pr( logFile, "Electrostatic energies will be calculated for all non-bonds between moving atoms.\n\n");
  868.         }
  869.         retval = sscanf( line, "%*s " FDFMT, &AD3_FE_coeff_estat );
  870.         if (retval == 1) {
  871.             if (outlev >= 0) {
  872.                 pr(logFile, "NOTE!  Internal electrostatics will NOT be scaled by the factor specified by this command,  %.4f -- the coefficient set by this command is ignored in AutoDock 4;\n", AD3_FE_coeff_estat);
  873.                 pr(logFile, "       the coefficient that will actually be used should be set in the parameter library file.\n");
  874.                 pr(logFile, "       The coefficient for the electrostatic energy term is %.4f", AD4.coeff_estat);
  875.                 if (parameter_library_found == 1) {
  876.                     pr( logFile, " as specified in parameter library \"%s\".\n", FN_parameter_library );
  877.                 } else {
  878.                     pr( logFile, ", the factory default value.\n");
  879.                 }
  880.             }
  881.         } else {
  882.             AD3_FE_coeff_estat = 1.0;
  883.         }
  884.  
  885.         (void) fflush(logFile);
  886.         break;
  887.  
  888. //______________________________________________________________________________
  889.  
  890.     case DPF_SEED:
  891.         /*
  892.         **  seed
  893.         **  Set the random-number gerator's seed value,
  894.         */
  895.         retval = (int)sscanf( line, "%*s %s %s", param[0], param[1]);
  896.         timeSeedIsSet[0] = 'F';
  897.         timeSeedIsSet[1] = 'F';
  898.         pr(logFile, "%d seed%c found.\n", retval, ((retval==1)? ' ' : 's'));
  899.         for (j=0; j<retval; j++) {
  900.             for (i=0; i<(int)strlen(param[j]); i++) {
  901.                 param[j][i] = (char)tolower( (int)param[j][i] );
  902.             }
  903.             pr(logFile, "argument \"%s\" found\n", param[j]);
  904.         }
  905.         if ((retval==2) || (retval==1)) {
  906.             for (i=0; i<retval ; i++ ) {
  907.                 if (equal(param[i], "tim", 3)) {
  908.                     timeSeedIsSet[i] = 'T';
  909.                     seed[i] = (FourByteLong)time( &time_seed );
  910.                     seed_random(seed[i]);
  911.                     pr(logFile,"Random number generator was seeded with the current time, value = %ld\n",seed[i]);
  912.                 } else if (equal(param[i], "pid", 3)) {
  913.                     timeSeedIsSet[i] = 'F';
  914.                     seed[i] = getpid();
  915.                     seed_random(seed[i]);
  916.                     pr(logFile,"Random number generator was seeded with the process ID, value   = %ld\n",seed[i]);
  917.                 } else {
  918.                     timeSeedIsSet[i] = 'F';
  919.                     seed[i] = atol(param[i]);
  920.                     seed_random(seed[i]);
  921.                     pr(logFile,"Random number generator was seeded with the user-specified value  %ld\n",seed[i]);
  922.                 }
  923.             }/*i*/
  924.             pr(logFile, "\n");
  925.             if (retval==2) {
  926.                 setall(seed[0], seed[1]);
  927.                 initgn(-1);  // Reinitializes the state of the current random number generator
  928.                 pr(logFile,"Portable random number generator was seeded with the user-specified values  %ld, %ld\n", seed[0], seed[1]);
  929.             }
  930.         } else {
  931.             pr(logFile, "Error encountered reading seeds!\n");
  932.         }
  933.  
  934.         (void) fflush(logFile);
  935.         break;
  936.  
  937. /*____________________________________________________________________________*/
  938.  
  939.     case DPF_LIGAND_TYPES:
  940.         /*
  941.          *  Read in the ligand atom type names, e.g.
  942.          *
  943.          *  ligand_types C HD OA P               # ligand atom type names
  944.          *
  945.          *  The order of the arguments is the index that will
  946.          *  be used for look up in the grid maps, "map_index".
  947.          */
  948.        
  949.         //  Use the function "parsetypes" to read in the atom types;
  950.         //  
  951.         //  The array "ligand_atom_type_ptrs" is returned, having been filled with pointers
  952.         //  to the beginning of each "atom type word" (not atom type characters);
  953.         //  In AutoDock 4, an atom type can be either 1 or 2 characters long.
  954.         //
  955.         //  Note: "atm_typ_str" (in AD3) served a similar role to "atom_type_name" (now used in AD4).
  956.         num_atom_types = parsetypes(line, ligand_atom_type_ptrs, MAX_ATOM_TYPES);
  957.  
  958.         B_found_ligand_types = TRUE;
  959.  
  960.         // This is not necessary if we increment num_maps one-at-a-time as we read each atom map in
  961.         // num_maps += num_atom_types;
  962.         info->num_atom_types = num_atom_types;
  963.  
  964.         for (i=0; i<num_atom_types; i++) {
  965.             strcpy(info->atom_type_name[i], ligand_atom_type_ptrs[i]);
  966. #ifdef DEBUG
  967.             (void) fprintf(logFile, "%d %s ->%s\n",i, ligand_atom_type_ptrs[i], info->atom_type_name[i]);
  968. #endif
  969.         }
  970.  
  971.         if (num_atom_types > 0) {
  972.             B_atom_types_found = TRUE;
  973.         } else {
  974.             prStr( error_message, "%s:  ERROR! No atom types have been found; we cannot continue without this information!\n\n", programname );
  975.             pr_2x( logFile, stderr, error_message );
  976.             prStr( error_message, "%s:  ERROR! Are you trying to use an AutoDock 3 DPF with AutoDock 4?\n\n", programname );
  977.             pr_2x( logFile, stderr, error_message );
  978.             exit(-1);
  979.         }
  980.  
  981.         if (debug > 0) {
  982.             for (i=0; i<num_atom_types; i++) {
  983.                 (void) fprintf(logFile, "info->atom_type_name[%d] = %s\n", i, info->atom_type_name[i] );
  984.             }
  985.         }
  986.  
  987.         // For all ligand atom types... set up the map_index
  988.         // "ligand_types"
  989.         for (i=0; i<num_atom_types; i++) {
  990.             foundParameter = apm_find(info->atom_type_name[i]);
  991.             if (foundParameter != NULL ) {
  992.                 // Not NULL means we have found this atom type's parameters.
  993.                 // Set the ParameterEntry's "map_index" member to the
  994.                 // 0-based index it had in the list of ligand types supplied in the DPF "types" line:
  995.                 foundParameter->map_index = i;
  996.                 parameterArray[i] = *(foundParameter);
  997.                 if (outlev > 0) {
  998.                     (void) fprintf( logFile, "Parameters found for ligand type \"%s\" (grid map index = %d, weighted well depth, epsilon = %6.4f)", foundParameter->autogrid_type, foundParameter->map_index, foundParameter->epsij );
  999.                     if (parameter_library_found == 1) {
  1000.                         pr( logFile, " in parameter library \"%s\".\n", FN_parameter_library );
  1001.                     } else {
  1002.                         pr( logFile, "\n");
  1003.                     }
  1004.                 }
  1005.             } else {
  1006.                 // We could not find this parameter -- return error here
  1007.                 prStr( error_message,"%s: ERROR:  Unknown ligand atom type \"%s\"; add parameters for it to the parameter library first!\n", programname, info->atom_type_name[i]);
  1008.                 pr_2x( logFile, stderr, error_message );
  1009.                 if (parameter_library_found == 1) {
  1010.                     prStr( error_message,"%s:         Edit the parameter library file \"%s\" and try again.\n", programname, FN_parameter_library );
  1011.                     pr_2x( logFile, stderr, error_message );
  1012.                 }
  1013.                 exit(-1);
  1014.             } // if / else apm_find
  1015.         } // for i
  1016.         pr( logFile, "\n\n");
  1017.  
  1018.         (void) fflush( logFile);
  1019.  
  1020.         // Calculate the internal energy table
  1021.  
  1022.         // loop over atom types, i, from 1 to number of atom types
  1023.         for (i=0; i<num_atom_types; i++) {
  1024.  
  1025.             //  Find internal energy parameters, i.e.  epsilon and r-equilibrium values...
  1026.             //  Lennard-Jones and Hydrogen Bond Potentials
  1027.  
  1028.             Ri = parameterArray[i].Rij;
  1029.             epsi = parameterArray[i].epsij;
  1030.             Ri_hb = parameterArray[i].Rij_hb;
  1031.             epsi_hb = parameterArray[i].epsij_hb;
  1032.             hbondi = parameterArray[i].hbond;
  1033.  
  1034.             // loop over atom types, j, from i to number of atom types
  1035.             for (j=i; j<num_atom_types; j++) {
  1036.                
  1037.                 //  Find internal energy parameters, i.e.  epsilon and r-equilibrium values...
  1038.                 //  Lennard-Jones and Hydrogen Bond Potentials
  1039.  
  1040.                 Rj = parameterArray[j].Rij;
  1041.                 epsj = parameterArray[j].epsij;
  1042.                 Rj_hb = parameterArray[j].Rij_hb;
  1043.                 epsj_hb = parameterArray[j].epsij_hb;
  1044.                 hbondj = parameterArray[j].hbond;
  1045.  
  1046.                 // we need to determine the correct xA and xB exponents
  1047.                 xA = 12; // for both LJ, 12-6 and HB, 12-10, xA is 12
  1048.                 xB =  6; // assume we have LJ, 12-6
  1049.  
  1050.                 if ( ((hbondi == DS) || (hbondi == D1)) && ((hbondj == AS) || (hbondj == A1) || (hbondj == A2)) ) {
  1051.                     // i is a donor and j is an acceptor.
  1052.                     // i is a hydrogen, j is a heteroatom
  1053.                     // we need to calculate the arithmetic mean of Ri_hb and Rj_hb  // not in this Universe...  :-(
  1054.                     //Rij = arithmetic_mean(Ri_hb, Rj_hb);// not in this Universe...  :-(
  1055.                     Rij = Rj_hb;
  1056.                     // we need to calculate the geometric mean of epsi_hb and epsj_hb  // not in this Universe...  :-(
  1057.                     //epsij = geometric_mean(epsi_hb, epsj_hb);// not in this Universe...  :-(
  1058.                     epsij = epsj_hb;
  1059.                     xB = 10;
  1060.                 } else if ( ((hbondi == AS) || (hbondi == A1) || (hbondi == A2)) && ((hbondj == DS) || (hbondj == D1))) {
  1061.                     // i is an acceptor and j is a donor.
  1062.                     // i is a heteroatom, j is a hydrogen
  1063.                     // we need to calculate the arithmetic mean of Ri_hb and Rj_hb// not in this Universe...  :-(
  1064.                     //Rij = arithmetic_mean(Ri_hb, Rj_hb);// not in this Universe...  :-(
  1065.                     Rij = Ri_hb;
  1066.                     // we need to calculate the geometric mean of epsi_hb and epsj_hb// not in this Universe...  :-(
  1067.                     //epsij = geometric_mean(epsi_hb, epsj_hb);// not in this Universe...  :-(
  1068.                     epsij = epsi_hb;
  1069.                     xB = 10;
  1070.                 } else {
  1071.                     // we need to calculate the arithmetic mean of Ri and Rj
  1072.                     Rij = arithmetic_mean(Ri, Rj);
  1073.                     // we need to calculate the geometric mean of epsi and epsj
  1074.                     epsij = geometric_mean(epsi, epsj);
  1075.                 }
  1076.  
  1077.                 /* Check that the Rij is reasonable */
  1078.                 if ((Rij < RIJ_MIN) || (Rij > RIJ_MAX)) {
  1079.                     (void) fprintf( logFile,
  1080.                     "WARNING: pairwise distance, Rij, %.2f, is not a very reasonable value for the equilibrium separation of two atoms! (%.2f Angstroms <= Rij <= %.2f Angstroms)\n\n", Rij, RIJ_MIN, RIJ_MAX);
  1081.                     (void) fprintf( logFile, "Perhaps you meant to use \"intnbp_coeffs\" instead of \"intnbp_r_eps\"?\n\n");
  1082.                     /* gmm commented out for dave goodsell, mutable atoms
  1083.                      * exit(-1); */
  1084.                 }
  1085.                 /* Check that the epsij is reasonable */
  1086.                 if ((epsij < EPSIJ_MIN) || (epsij > EPSIJ_MAX)) {
  1087.                     (void) fprintf( logFile,
  1088.                     "WARNING: well-depth, epsilon_ij, %.2f, is not a very reasonable value for the equilibrium potential energy of two atoms! (%.2f kcal/mol <= epsilon_ij <= %.2f kcal/mol)\n\n", epsij, EPSIJ_MIN, EPSIJ_MAX);
  1089.                     (void) fprintf( logFile, "Perhaps you meant to use \"intnbp_coeffs\" instead of \"intnbp_r_eps\"?\n\n");
  1090.                     /* gmm commented out for dave goodsell, mutable atoms
  1091.                      * exit(-1); */
  1092.                 }
  1093.                 /* Defend against division by zero... */
  1094.                 if (xA != xB) {
  1095.                     cA = (tmpconst = epsij / (Real)(xA - xB)) * pow( (double)Rij, (double)xA ) * (Real)xB;
  1096.                     cB = tmpconst * pow( (double)Rij, (double)xB ) * (Real)xA;
  1097.                     pr(logFile, "\nCalculating internal non-bonded interaction energies for docking calculation;\n");
  1098.                     intnbtable( &B_havenbp, a1, a2, info, cA, cB, xA, xB, AD4.coeff_desolv, sigma, ad_energy_tables, BOUND_CALCULATION );
  1099.                     pr(logFile, "\nCalculating internal non-bonded interaction energies for unbound conformation calculation;\n");
  1100.                     intnbtable( &B_havenbp, a1, a2, info, cA_unbound, cB_unbound, xA_unbound, xB_unbound, AD4.coeff_desolv, sigma, unbound_energy_tables, UNBOUND_CALCULATION );
  1101.                     // Increment the atom type numbers, a1 and a2, for the internal non-bond table
  1102.                     a2++;
  1103.                     if (a2 >= info->num_atom_types) {
  1104.                         a1++;
  1105.                         a2 = a1;
  1106.                     }
  1107.                    
  1108.                 } else {
  1109.                     pr(logFile,"WARNING: Exponents must be different, to avoid division by zero!\n\tAborting...\n");
  1110.                     exit(-1);
  1111.                 }
  1112.                 (void) fflush(logFile);
  1113.  
  1114.             } // for j
  1115.         } // for i
  1116.         break;
  1117.  
  1118. //______________________________________________________________________________
  1119.  
  1120.     case DPF_FLD:
  1121.         /*
  1122.         ** fld
  1123.         ** GRID_DATA_FILE
  1124.         ** Read the (AVS-format) grid data file, .fld
  1125.         */
  1126.         // TO DO: add outlev
  1127.         readfield( info, line, jobStart, tms_jobStart );
  1128.  
  1129.         /*
  1130.         // Dynamically allocate memory for the maps
  1131.         map = NewGridMapSet(info);
  1132.            
  1133.         if (map == NULL) {
  1134.             prStr(error_message, "%s:  Sorry, there is not enough memory to store the grid maps.  Please use smaller maps and/or fewer atom types.\n", programname);
  1135.             stop(error_message);
  1136.             exit(1);
  1137.         }
  1138.         // Initialise the maps
  1139.         for (i=0; i<num_map_values; i++) {
  1140.             map[i] = 0.0L;
  1141.         }
  1142.         */
  1143.         (void) fflush(logFile);
  1144.         break;
  1145.  
  1146. //______________________________________________________________________________
  1147.  
  1148.     case DPF_MAP:
  1149.         /*
  1150.         ** map
  1151.         ** ATOMIC AFFINITY, ELECTROSTATIC POTENTIAL OR DESOLVATION ENERGY GRID MAP
  1152.         ** Read in active site grid map...
  1153.         */
  1154.         B_charMap = FALSE;
  1155.         if (B_atom_types_found == TRUE) {
  1156.             // Read in the AutoGrid atomic affinity map
  1157.             // map_index could be incremented here if we had the atom_type stored in each map...
  1158.             map_stats = readmap( line, outlev, jobStart, tms_jobStart, B_charMap, &B_havemap, num_maps, info, map, 'a' );
  1159.             pr(logFile, "Min= %.3lf Mean= %.3lf Max= %.3lf\n\n",
  1160.                     map_stats.minimum, map_stats.mean, map_stats.maximum);
  1161.             num_maps++;
  1162.  
  1163.         } else {
  1164.             prStr( error_message, "%s:  ERROR! No atom types have been found; we cannot continue without this information!\n\n", programname );
  1165.             pr_2x( logFile, stderr, error_message );
  1166.             prStr( error_message, "%s:  ERROR! Are you trying to use an AutoDock 3 DPF with AutoDock 4?\n\n", programname );
  1167.             pr_2x( logFile, stderr, error_message );
  1168.             exit(-1);
  1169.         }
  1170.         (void) fflush(logFile);
  1171.         break;
  1172.  
  1173. //______________________________________________________________________________
  1174.  
  1175.     case DPF_ELECMAP:
  1176.         /*
  1177.          *  elecmap file.e.map
  1178.          */
  1179.         map_stats = readmap( line, outlev, jobStart, tms_jobStart, B_charMap, &B_havemap, num_maps, info, map, 'e' );
  1180.         pr(logFile, "Min= %.3lf Mean= %.3lf Max= %.3lf\n\n",
  1181.                 map_stats.minimum, map_stats.mean, map_stats.maximum);
  1182.         ElecMap = num_maps;
  1183.         B_found_elecmap = TRUE;
  1184.         num_maps++;
  1185.         (void) fflush(logFile);
  1186.         break;
  1187.  
  1188. //______________________________________________________________________________
  1189.  
  1190.     case DPF_DESOLVMAP:
  1191.         /*
  1192.          *  desolvmap file.d.map
  1193.          */
  1194.         map_stats = readmap( line, outlev, jobStart, tms_jobStart, B_charMap, &B_havemap, num_maps, info, map, 'd' );
  1195.         pr(logFile, "Min= %.3lf Mean= %.3lf Max= %.3lf\n\n",
  1196.                 map_stats.minimum, map_stats.mean, map_stats.maximum);
  1197.         DesolvMap = num_maps;
  1198.         B_found_desolvmap = TRUE;
  1199.         num_maps++;
  1200.         (void) fflush(logFile);
  1201.         break;
  1202.  
  1203. //______________________________________________________________________________
  1204.  
  1205.     case DPF_CHARMAP:
  1206.         /*
  1207.         ** charmap
  1208.         ** ATOMIC AFFINITY, ELECTROSTATIC POTENTIAL OR DESOLVATION ENERGY GRID MAP
  1209.         ** Read in active site grid map...
  1210.         */
  1211.         B_charMap = TRUE;
  1212.         if (B_atom_types_found == TRUE) {
  1213.             // map_index could be incremented here if we had the atom_type stored in each map...
  1214.             map_stats = readmap( line, outlev, jobStart, tms_jobStart, B_charMap, &B_havemap, num_maps, info, map, 'c' );
  1215.             pr(logFile, "Min= %.3lf Mean= %.3lf Max= %.3lf\n\n",
  1216.                     map_stats.minimum, map_stats.mean, map_stats.maximum);
  1217.             num_maps++;
  1218.         } else {
  1219.             prStr( error_message, "%s:  ERROR! No atom types have been found; we cannot continue without this information!\n\n", programname );
  1220.             pr_2x( logFile, stderr, error_message );
  1221.             prStr( error_message, "%s:  ERROR! Are you trying to use an AutoDock 3 DPF with AutoDock 4?\n\n", programname );
  1222.             pr_2x( logFile, stderr, error_message );
  1223.             exit(-1);
  1224.         }
  1225.         (void) fflush(logFile);
  1226.         break;
  1227.  
  1228. //______________________________________________________________________________
  1229.  
  1230.     case DPF_MOVE:
  1231.         /*
  1232.         ** move ligand_file.pdbqt
  1233.         ** Specify the movable ligand,
  1234.         */
  1235.         //
  1236.         // Initialisations that must be done before reading in a new ligand...
  1237.         //
  1238.         nconf = 0;
  1239.         for (k = 0; k < MAX_RUNS; k++) {
  1240.             for (i = 0; i  < MAX_TORS;  i++ ) {
  1241.                 sHist[k].tor[i] = 0.0;
  1242.             }
  1243.         }
  1244.         for (j = 0; j < MAX_RUNS; j++) {
  1245.             // isort[j] = j;
  1246.             econf[j] = 0.0;
  1247.         }
  1248.         for (j = 0;  j < MAX_ATOMS;  j++ ) {
  1249.             type[j] = 0;
  1250.             ignore_inter[j] = 0;
  1251.         }
  1252.         for (i = 0; i  < MAX_TORS;  i++ ) {
  1253.             for (j = 0;  j < MAX_ATOMS;  j++ ) {
  1254.                 tlist[i][j] = 0;
  1255.             }
  1256.         }
  1257.         for (i = 0; i  < MAX_TORS;  i++ ) {
  1258.             B_isTorConstrained[i] = 0;
  1259.             US_torProfile[i][0] = 0;
  1260.             N_con[i] = 0;
  1261.         }
  1262.         initialiseState( &sInit );
  1263.         initialiseState( &(ligand.S) );
  1264.         initialiseQuat( &q_reorient );
  1265.         B_constrain_dist = B_haveCharges = FALSE;
  1266.         ntor1 = ntor = atomC1 = atomC2 = 0;
  1267.         sqlower = squpper = 0.0;
  1268.         strcpy( FN_pop_file, "");  // means don't print pop_file
  1269.         //
  1270.         // end of initialization
  1271.         //
  1272.  
  1273.         // this is the DPF_MOVE section...
  1274.         B_found_move_keyword = TRUE;
  1275.  
  1276.         print_1_4_message(logFile, B_include_1_4_interactions, scale_1_4);
  1277.  
  1278.         natom=0;
  1279.         ligand = readPDBQT( line,
  1280.                             num_atom_types,
  1281.                             &natom,
  1282.                             crdpdb, crdreo, charge, &B_haveCharges,
  1283.                             type, bond_index,
  1284.                             pdbaname, FN_ligand, FN_flexres, B_have_flexible_residues, atomstuff, Htype,
  1285.                             &B_constrain_dist, &atomC1, &atomC2,
  1286.                             &sqlower, &squpper,
  1287.                             &ntor1, &ntor, &ntor_ligand,
  1288.                             tlist, vt,
  1289.                             &Nnb, nonbondlist,
  1290.                             jobStart, tms_jobStart, hostnm, &ntorsdof, outlev,
  1291.                             ignore_inter,
  1292.                             B_include_1_4_interactions,
  1293.                             atoms, PDBQT_record );
  1294.  
  1295.         // pre-calculate some values we will need later in computing the desolvation energy
  1296.         //
  1297.         for (i=0; i<natom; i++) {
  1298.             abs_charge[i] = fabs(charge[i]);
  1299.             qsp_abs_charge[i] = qsolpar * abs_charge[i];
  1300.         }
  1301.         pr(logFile, "Number of atoms in ligand:  %d\n\n", true_ligand_atoms);
  1302.  
  1303.         pr(logFile, "Number of vibrational degrees of freedom of ligand:  %d\n\n\n", (3 * true_ligand_atoms) - 6 );
  1304.         pr( logFile, "Number of torsional degrees of freedom = %d\n", ntorsdof);
  1305.  
  1306.         torsFreeEnergy = (Real)ntorsdof * AD4.coeff_tors;
  1307.  
  1308.         pr( logFile, "Estimated loss of torsional free energy upon binding = %+.4f kcal/mol\n\n", torsFreeEnergy);
  1309.  
  1310.         for (i=0;i<natom;i++) {
  1311.             if (ignore_inter[i] == 1) {
  1312.                 pr(logFile, "Special Boundary Conditions:\n");
  1313.                 pr(logFile, "____________________________\n\n");
  1314.                 pr(logFile, "AutoDock will ignore the following atoms in the input PDBQT file \nin intermolecular energy calculations:\n");
  1315.                 pr(logFile, "\n(This is because these residue atoms are at the boundary between \nflexible and rigid, and since they cannot move \nthey will not affect the total energy.)\n\n");
  1316.                 break;
  1317.             }
  1318.         }
  1319.         for (i=0;i<natom;i++) {
  1320.             if (ignore_inter[i] == 1) {
  1321.                 pr(logFile, "Atom number %d:  %s\n", i+1, atomstuff[i] );
  1322.             }
  1323.         }
  1324.         pr(logFile, "\n");
  1325.  
  1326.         if (!B_haveCharges) {
  1327.             pr( logFile, "%s: WARNING! No partial atomic charges have been supplied yet.\n\n",programname);
  1328.         } else {
  1329.             if (Nnb > 0) {
  1330.             pr(logFile,"Calculating the product of the partial atomic charges, q1*q2, for all %d non-bonded pairs...\n\n", Nnb);
  1331.             pr(logFile," -- Scaled q1*q2 means multiplied by both  %.1lf (for conversion later on to kcal/mol)\n", (double)ELECSCALE);
  1332.             pr(logFile,"    and by the AD4 FF electrostatics coefficient, %.4lf\n\n", (double)AD4.coeff_estat);
  1333.             if (outlev >= 0) {
  1334.                 pr(logFile,"Non-bonded                           Scaled\n");
  1335.                 pr(logFile,"   Pair     Atom1-Atom2    q1*q2      q1*q2\n");
  1336.                 pr(logFile,"__________  ___________  _________  _________\n");
  1337.                 for (i = 0;  i < Nnb;  i++) {
  1338.                     atm1 = nonbondlist[i].a1;
  1339.                     atm2 = nonbondlist[i].a2;
  1340.                     int t1 = nonbondlist[i].t1;
  1341.                     int t2 = nonbondlist[i].t2;
  1342.                     nonbondlist[i].desolv =
  1343.                            ( parameterArray[t2].vol * (parameterArray[t1].solpar + qsp_abs_charge[atm1])
  1344.                            + parameterArray[t1].vol * (parameterArray[t2].solpar + qsp_abs_charge[atm2]) );
  1345.                     nonbondlist[i].q1q2 = charge[atm1] * charge[atm2];
  1346.                     pr(logFile,"   %4d     %5d-%-5d    %5.2f",i+1,atm1+1,atm2+1,nonbondlist[i].q1q2);
  1347.                     nonbondlist[i].q1q2 *= ELECSCALE * AD4.coeff_estat;
  1348.                     pr(logFile,"     %6.2f\n",nonbondlist[i].q1q2);
  1349.                 }
  1350.                 pr(logFile,"\n");
  1351.                 } // if outlev > 0
  1352.             } // if NNb > 0
  1353.         } // else
  1354.  
  1355.         sInit.ntor = ligand.S.ntor;
  1356.         ++nmol;
  1357.         ++nlig;
  1358.  
  1359.         (void) fflush(logFile);
  1360.         break;
  1361.  
  1362. /*____________________________________________________________________________*/
  1363.  
  1364.     case DPF_FLEXRES:
  1365.         /*
  1366.          * flexible_residues file.pdbqt
  1367.          */
  1368.         (void) fflush(logFile);
  1369.         pr(logFile, "\nThe flexible residues will be read in from \"%s\".\n", FN_flexres);
  1370.         break;
  1371.  
  1372.  
  1373. #ifdef USING_COLINY
  1374. /*____________________________________________________________________________*/
  1375.  
  1376.     case DPF_COLINY:
  1377.     {
  1378.         //ostdiostream fstr(logFile);
  1379.         //ostdiostream fstr(logFile->_file);
  1380.         //CommonIO::set_streams(&fstr,&fstr,&cin);
  1381.  
  1382.         struct tms tms_colinyStart;
  1383.         struct tms tms_colinyEnd;
  1384.  
  1385.         Clock  colinyStart;
  1386.         Clock  colinyEnd;
  1387.  
  1388.         int coliny_seed;
  1389.         char algname[64];
  1390.         char nruns_str[64];
  1391.         (void) sscanf(line, "%*s %s %d", algname, &nruns);
  1392.         (void) sscanf(line, "%*s %s %s", algname, nruns_str);
  1393.  
  1394.         if (strcmp(algname,"help")==0) {
  1395.             std::vector<double> initvec;
  1396.             coliny_init(algname, "");
  1397.             prStr(error_message, "%s:  ERROR:  no optimizer type specified.", programname);
  1398.             stop(error_message);
  1399.             exit(-1);
  1400.         }
  1401.         else if (strcmp(nruns_str,"help")==0) {
  1402.             std::vector<double> initvec;
  1403.             coliny_init(algname, nruns_str);
  1404.             prStr(error_message, "%s:  ERROR:  no optimizer type specified.", programname);
  1405.             stop(error_message);
  1406.             exit(-1);
  1407.         }
  1408.  
  1409.         if (!command_mode) {
  1410.             if (nruns>MAX_RUNS) {
  1411.                 prStr(error_message, "%s:  ERROR:  %d runs requested, but only dimensioned for %d.\nChange \"MAX_RUNS\" in \"constants.h\".", programname, nruns, MAX_RUNS);
  1412.                 stop(error_message);
  1413.                 exit(-1);
  1414.             } else if (!B_found_elecmap) {
  1415.                 prStr(error_message, "%s:  ERROR:  no \"elecmap\" command has been specified!\n", programname);
  1416.                 stop(error_message);
  1417.                 exit(-1);
  1418.             } else if (!B_found_desolvmap) {
  1419.                 prStr(error_message, "%s:  ERROR:  no \"desolvmap\" command has been specified!\n", programname);
  1420.                 stop(error_message);
  1421.                 exit(-1);
  1422.             }
  1423.  
  1424.             evaluate.setup( crd, charge, abs_charge, qsp_abs_charge, type, natom,
  1425.                             map, elec, emap, nonbondlist, ad_energy_tables,
  1426.                             Nnb, B_calcIntElec, B_isGaussTorCon, B_isTorConstrained, B_ShowTorE,
  1427.                             US_TorE, US_torProfile,
  1428.                             vt, tlist,
  1429.                             crdpdb, crdreo, sInit, ligand, ignore_inter, B_include_1_4_interactions, scale_1_4,
  1430.                             parameterArray, unbound_internal_FE, info,
  1431.                             B_use_non_bond_cutoff, B_have_flexible_residues);
  1432.  
  1433.             evaluate.compute_intermol_energy(TRUE);
  1434.  
  1435.             char domain[1024];
  1436.             // NOTE: Coliny enforces the bound constraints, but since the
  1437.             // torsion angles are periodic, we simply prevent the optimizer
  1438.             // from going too far.
  1439.             if (sInit.ntor > 0) {
  1440.                 sprintf(domain,"[%f,%f] [%f,%f] [%f,%f] [-1000.0,1000.0]^3 [-3.1416,3.1416] [-3.1416,3.1416]^%d",(double)info->lo[X], (double)info->hi[X], (double)info->lo[Y], (double)info->hi[Y], (double)info->lo[Z], (double)info->hi[Z], sInit.ntor);
  1441.                 //sprintf(domain,"[%f,%f] [%f,%f] [%f,%f] [-1.0,1.1]^3 [-6.2832,12.5664] [-6.2832,12.5664]^%d",(double)info->lo[X], (double)info->hi[X], (double)info->lo[Y], (double)info->hi[Y], (double)info->lo[Z], (double)info->hi[Z], sInit.ntor);
  1442.             } else {
  1443.                 sprintf(domain,"[%f,%f] [%f,%f] [%f,%f] [-1000.0,1000.0]^3 [-3.1416,3.1416]",(double)info->lo[X], (double)info->hi[X], (double)info->lo[Y], (double)info->hi[Y], (double)info->lo[Z], (double)info->hi[Z]);
  1444.             }
  1445.             pr(logFile, "Number of Coliny %s dockings = %d run%c\n", algname, nruns, (nruns>1)?'s':' ');
  1446.             pr(logFile, "Search Domain: %s\n", domain);
  1447.  
  1448.             //
  1449.             // COLINY-SPECIFIC LOGIC - BEGIN
  1450.             //
  1451.  
  1452.             try {
  1453.  
  1454.                 std::vector<double> initvec, finalpt;
  1455.                 // set up initial point
  1456.                 initvec.resize(7+sInit.ntor);
  1457.                 initvec[0] = sInit.T.x;
  1458.                 initvec[1] = sInit.T.y;
  1459.                 initvec[2] = sInit.T.z;
  1460.                 /*
  1461.                  * axis-angle (nx,ny,nz,ang) suffers from bias
  1462.                 initvec[3] = sInit.Q.nx;
  1463.                 initvec[4] = sInit.Q.ny;
  1464.                 initvec[5] = sInit.Q.nz;
  1465.                 initvec[6] = DegreesToRadians( sInit.Q.ang );
  1466.                 */
  1467.                 sInit.Q = convertRotToQuat( sInit.Q );
  1468.                 initvec[3] = sInit.Q.x;
  1469.                 initvec[4] = sInit.Q.y;
  1470.                 initvec[5] = sInit.Q.z;
  1471.                 initvec[6] = sInit.Q.w;
  1472.                 for (j=0; j < sInit.ntor ; j++) {
  1473.                   initvec[j+7] = DegreesToRadians(sInit.tor[j]);
  1474.                 }
  1475.                 coliny_init(algname, domain);
  1476.  
  1477.                 for (j=0; j<nruns; j++) {
  1478.                   fprintf( logFile, "\n\n\tBEGINNING Coliny %s DOCKING\n",algname);
  1479.                   pr(logFile, "\nDoing %s run:  %d/%d.\n", algname, j+1, nruns);
  1480.  
  1481.                   //coliny uses a single seed
  1482.                   coliny_seed = seed[0]+seed[1]+j;
  1483.                   pr(logFile, "Seed: %d [%ld+%ld+%d]\n", coliny_seed, seed[0], seed[1], j);
  1484.                   //pr(logFile, "Seeds:  %ld %ld\n", seed[0], seed[1]);
  1485.                   (void) fflush(logFile);
  1486.  
  1487.                   colinyStart = times(&tms_colinyStart);
  1488.  
  1489.                   finalpt.resize( initvec.size() );
  1490.                   int neval, niters;
  1491.                   coliny_minimize( coliny_seed, initvec, finalpt, neval, niters );
  1492.                   //fstr.flush();
  1493.  
  1494.                   make_state_from_rep( (double *)&(finalpt[0]), int(finalpt.size()), &sHist[nconf]);
  1495.  
  1496.                   pr(logFile, "\nTotal Num Evals: %d\n", neval);
  1497.                   printState(logFile, sHist[nconf], 2);
  1498.  
  1499.                   colinyEnd = times(&tms_colinyEnd);
  1500.                   pr(logFile, "Time taken for this %s run:\n", algname);
  1501.                   timesyshms(colinyEnd-colinyStart, &tms_colinyStart, &tms_colinyEnd);
  1502.                   pr(logFile, "\n");
  1503.  
  1504.                   pr(logFile, "Total number of Energy Evaluations: %d\n", (int)evaluate.evals() );
  1505.                   //pr(logFile, "Total number of Iterations:        %d\n", (int)niters);
  1506.  
  1507.                   pr(logFile, "\nFinal docked state:\n");
  1508.                   pr( logFile, UnderLine );
  1509.                   pr( logFile, "\n\n\tFINAL Coliny %s DOCKED STATE\n",algname );
  1510.                   pr( logFile,     "\t____________________________________\n\n\n" );
  1511.                   (void) fflush(logFile);
  1512.  
  1513.                   writePDBQT( j, seed, FN_ligand, dock_param_fn, lig_center,
  1514.                               sHist[nconf], ntor, &eintra, &einter, natom, atomstuff,
  1515.                               crd, emap, elec,
  1516.                               charge, abs_charge, qsp_abs_charge,
  1517.                               ligand_is_inhibitor,
  1518.                               torsFreeEnergy,
  1519.                               vt, tlist, crdpdb, nonbondlist,
  1520.                               ad_energy_tables,
  1521.                               type, Nnb, B_calcIntElec,
  1522.                               map,
  1523.                               outlev,
  1524.                               ignore_inter,
  1525.                               B_include_1_4_interactions, scale_1_4, parameterArray, unbound_internal_FE,
  1526.                               info, DOCKED, PDBQT_record, B_use_non_bond_cutoff, B_have_flexible_residues);
  1527.  
  1528.                   econf[nconf] = eintra + einter + torsFreeEnergy - unbound_internal_FE;
  1529.                   evaluate.reset();
  1530.                  
  1531.                   ++nconf;
  1532.  
  1533.                 } // Next run
  1534.                 if(write_stateFile){
  1535.                   fprintf(stateFile,"\t</runs>\n");
  1536.                   (void) fflush(stateFile);
  1537.                 }
  1538.                 (void) fflush(logFile);
  1539.             }
  1540.             catch (std::exception& err) {
  1541.               (void)fprintf(logFile, "Caught Exception: %s\n", err.what());
  1542.               exit(1);
  1543.             }
  1544.  
  1545.         } else {
  1546.             (void)fprintf(logFile, "NOTE: Command mode has been set, so optimization cannot be performed.\n\n");
  1547.         }
  1548.     }
  1549.     break;
  1550. #endif
  1551.  
  1552.  
  1553. //______________________________________________________________________________
  1554.  
  1555.     case DPF_ABOUT:
  1556.         /*
  1557.         **  about
  1558.         **  Rotation center for current ligand,
  1559.         */
  1560.         (void) sscanf( line, "%*s " FDFMT3, &lig_center[X], &lig_center[Y], &lig_center[Z]);
  1561.         pr( logFile, "Small molecule center of rotation =\t" );
  1562.         pr( logFile, "(%+.3f, %+.3f, %+.3f)\n\n", lig_center[X], lig_center[Y], lig_center[Z]);
  1563.         /*
  1564.         **  Center the ligand,
  1565.         */
  1566.         if ( nmol == 0 ) {
  1567.             pr( logFile, "Must specify a ligand PDBQT file, using the \"move\" command.\n");
  1568.         } else {
  1569.             if (outlev >= 0) {
  1570.                 pr( logFile, "Translating small molecule by:\t" );
  1571.                 pr( logFile, "(%+.3f, %+.3f, %+.3f)\n\n", -lig_center[X], -lig_center[Y], -lig_center[Z]);
  1572.             }
  1573.             /*
  1574.             **  Zero-out on central point...
  1575.             */
  1576.             maxrad = -1.0;
  1577.             for ( i=0; i<true_ligand_atoms; i++ ) { /*new, gmm, 6-23-1998*/
  1578.                 r2sum=0.0;
  1579.                 for (xyz = 0;  xyz < SPACE;  xyz++) {
  1580.                     c = crd[i][xyz] = (crdpdb[i][xyz] -= lig_center[xyz]);
  1581.                     r2sum += c*c;
  1582.                 } /* xyz */
  1583.                 maxrad = max(maxrad,sqrt(r2sum));
  1584.             } /* i */
  1585.             if (outlev >= 0) {
  1586.                 pr( logFile, "Furthest ligand atom from \"about\" center is %.3f Angstroms (maxrad).\n\n",maxrad);
  1587.             }
  1588.         }
  1589.         (void) fflush(logFile);
  1590.         break;
  1591.  
  1592. /*____________________________________________________________________________*/
  1593.  
  1594.     case DPF_REORIENT:
  1595.         /*
  1596.          *  reorient random
  1597.          *      # applies a random rotation to the input ligand
  1598.          * -OR-
  1599.          *  reorient standard
  1600.          *      # moves the ligand such that
  1601.          *      # the first three atoms lie parallel to the xy-plane, and
  1602.          *      # the first two atoms lie parallel to the x-axis
  1603.          * -OR-
  1604.          *  reorient <axis-x> <axis-y> <axis-z> <angle>
  1605.          *      # applies the specified rotation to the input ligand
  1606.          */
  1607.         (void) sscanf( line, "%*s %s", param[0] );
  1608.         { // Parse the reorient command
  1609.             for (i=0; i<6; i++) {
  1610.                 param[0][i] = (char)tolower( (int)param[0][i] );
  1611.             }
  1612.             if (equal(param[0],"random",6)) {
  1613.                 // reorient random
  1614.                 B_reorient_random = TRUE; // create a new random orientation before docking
  1615.  
  1616.                 create_random_orientation( &q_reorient );
  1617.  
  1618.             } else if (equal(param[0],"standard",8)) {
  1619.                 { // reorient standard
  1620.                 B_reorient_random = FALSE; // do not create a new random orientation before docking
  1621.  
  1622.                 if (true_ligand_atoms >= 3 ) {
  1623.                     // Move the ligand such that
  1624.                     // the first three atoms lie parallel to the xy-plane, and
  1625.                     // the first two atoms lie parallel to the x-axis
  1626.                     Vector vec_01,     // vector between ligand atoms 0 and 1
  1627.                            vec_12,     // vector between ligand atoms 1 and 2
  1628.                            vec_normal, // vector perpendicular to plane of vec_01 and vec_12
  1629.                            vec_x_axis, // vector along the X-axis
  1630.                            vec_z_axis, // vector along the Z-axis
  1631.                            vec_reorient_axis; // vector describing the axis about which to reorient
  1632.                     // Set the X and Z axes:
  1633.                     vec_x_axis[X] = 1.;
  1634.                     vec_x_axis[Y] = 0.;
  1635.                     vec_x_axis[Z] = 0.;
  1636.                     vec_z_axis[X] = 0.;
  1637.                     vec_z_axis[Y] = 0.;
  1638.                     vec_z_axis[Z] = 1.;
  1639.                     for (xyz = 0;  xyz < SPACE;  xyz++) {
  1640.                         vec_01[xyz] = (double)( crdpdb[1][xyz] - crdpdb[0][xyz] );
  1641.                         vec_12[xyz] = (double)( crdpdb[2][xyz] - crdpdb[1][xyz] );
  1642.                     }
  1643.                     // Compute the normal to vec_01 and vec_12
  1644.                     Cross_product( vec_normal, vec_01, vec_12 );
  1645.                     Print_vector( logFile, "vec_01", vec_01 );
  1646.                     Print_vector( logFile, "vec_12", vec_12 );
  1647.                     Print_vector( logFile, "vec_normal", vec_normal );
  1648.                     Print_vector( logFile, "vec_z_axis", vec_z_axis );
  1649.                     // Compute the angle between vec_01 and vec_12
  1650.                     double angle_012 = 0.;
  1651.                     angle_012 = Angle_between( vec_01, vec_12 );
  1652.                     pr( logFile, "Angle between vectors 01 and 12 = %.2f degrees\n", RadiansToDegrees( angle_012 ) );
  1653.                     if ( ( fabs(angle_012) < APPROX_ZERO ) || ( ( fabs(angle_012) > (PI - APPROX_ZERO) ) && ( fabs(angle_012) < (PI + APPROX_ZERO) ) ) ) {
  1654.                         // angle is too small or "too linear" to align the molecule into the xy-plane
  1655.                         pr( logFile, "%s:  WARNING!  The angle between the first three atoms is not suitable (%6.3f degrees) to align them with the xy-plane.\n", programname, RadiansToDegrees( angle_012 ) );
  1656.                     } else {
  1657.                         // Calculate angle between vec_normal and the z-axis
  1658.                         double angle_n1z = 0.;  // Angle between vec_normal and the z-axis
  1659.                         angle_n1z = Angle_between( vec_normal, vec_z_axis );
  1660.                         pr( logFile, "Angle between vec_normal and vec_z_axis = %.2f degrees\n", RadiansToDegrees( angle_n1z ) );
  1661.                         //
  1662.                         // We need to rotate the molecule about the normal to vec_normal and vec_z_axis
  1663.                         Cross_product( vec_reorient_axis, vec_normal, vec_z_axis );
  1664.                         //
  1665.                         // Set the rotation axis for reorientation
  1666.                         q_reorient.nx = vec_reorient_axis[X];
  1667.                         q_reorient.ny = vec_reorient_axis[Y];
  1668.                         q_reorient.nz = vec_reorient_axis[Z];
  1669.                         //
  1670.                         // Normalise the vector defining the axis of rotation:
  1671.                         q_reorient = normRot( q_reorient );
  1672.                         //
  1673.                         // Set the angle for reorientation of the first 3 atoms
  1674.                         // into the xy-plane
  1675.                         q_reorient.ang = -angle_n1z;
  1676.                         //
  1677.                         // Convert the rotation-about-axis components (nx,ny,nz,ang)
  1678.                         // to a rotation-quaternion (x,y,z,w):
  1679.                         q_reorient = convertRotToQuat( q_reorient );
  1680.  
  1681.                         // Rotate ligand into the xy-plane...
  1682.                         // qtransform( origin, q_reorient, crdreo, true_ligand_atoms );
  1683.                         qtransform( origin, q_reorient, crdpdb, true_ligand_atoms );
  1684.  
  1685.                         // Compute the updated vec_01, the vector between atom 0 and atom 1,
  1686.                         // since the preceding "qtransform" changed the coordinates.
  1687.                         for (xyz = 0;  xyz < SPACE;  xyz++) {
  1688.                             // vec_01[xyz] = (double)( crdreo[1][xyz] - crdreo[0][xyz] );
  1689.                             vec_01[xyz] = (double)( crdpdb[1][xyz] - crdpdb[0][xyz] );
  1690.                         }
  1691.                         //
  1692.                         // Compute the angle between vec_01 and the x-axis:
  1693.                         double angle_01x = 0.;
  1694.                         angle_01x = Angle_between( vec_01, vec_x_axis );
  1695.                         //
  1696.                         pr( logFile, "Angle between vector 01 and the x-axis = %.2f degrees\n", RadiansToDegrees( angle_01x ) );
  1697.                         //
  1698.                         // The rotation axis to rotate the first two atoms, 0 and 1,
  1699.                         // to be parallel to the x-axis, will be
  1700.                         // perpendicular to the xy-plane, i.e. the z-axis,
  1701.                         // since the molecule's first 3 atoms are now in the xy-plane.
  1702.                         q_reorient.nx = vec_z_axis[X];
  1703.                         q_reorient.ny = vec_z_axis[Y];
  1704.                         q_reorient.nz = vec_z_axis[Z];
  1705.                         //
  1706.                         // Set the rotation angle:
  1707.                         q_reorient.ang = angle_01x;
  1708.                         //
  1709.                         // Build the quaternion from the axis-angle rotation values:
  1710.                         q_reorient = convertRotToQuat( q_reorient );
  1711.                     } // angle_012 is appropriate to align into xy-plane
  1712.  
  1713.                 } else {
  1714.                     prStr( error_message, "%s: ERROR! Insufficient atoms in the ligand.  There must be at least three atoms in the ligand to use this command.\n", programname );
  1715.                     stop( error_message );
  1716.                     exit( -1 );
  1717.                 }
  1718.                 } // reorient standard
  1719.             } else {
  1720.                 { // reorient <nx> <ny> <nz> <angle>
  1721.                     B_reorient_random = FALSE; // do not create a new random orientation before docking
  1722.  
  1723.                     // Read the specified initial orientation for the ligand
  1724.                     retval = (int)sscanf( line,"%*s %lf %lf %lf %lf", &(q_reorient.nx), &(q_reorient.ny), &(q_reorient.nz), &(q_reorient.ang) );
  1725.                     if ( retval == 4 ) {
  1726.                         // Normalise the vector defining the axis of rotation:
  1727.                         q_reorient = normRot( q_reorient );
  1728.                         // Make sure angle is in radians, and ranges from -PI to PI
  1729.                         q_reorient.ang = DegreesToRadians( q_reorient.ang ); // convert from degrees to radians
  1730.                         q_reorient.ang = ModRad( q_reorient.ang ); // wrap to range (0, 2*PI) using modulo 2*PI
  1731.                         q_reorient.ang = WrpRad( q_reorient.ang ); // wrap to range (-PI, PI)
  1732.                         pr( logFile, "After normalising the vector, and converting the angle to radians, the axis-angle rotation becomes ((%.3f, %.3f, %.3f), %.2f radians)\n",
  1733.                                 q_reorient.nx, q_reorient.ny, q_reorient.ny, q_reorient.ang);
  1734.                         // Convert the rotation-about-axis components (nx,ny,nz,ang)
  1735.                         // to a rotation-quaternion (x,y,z,w):
  1736.                         q_reorient = convertRotToQuat( q_reorient );
  1737.                     } else {
  1738.                         prStr( error_message, "%s: ERROR! Please specify the vector and rotation angle using four real numbers.\n", programname );
  1739.                         stop( error_message );
  1740.                         exit( -1 );
  1741.                     }
  1742.                 } // reorient <nx> <ny> <nz> <angle>
  1743.             } // endif
  1744.         } // end parsing reorient command line
  1745.  
  1746.         // reorient( logFile, true_ligand_atoms, atomstuff, crdreo, charge, type,
  1747.         reorient( logFile, true_ligand_atoms, atomstuff, crdpdb, charge, type,
  1748.                   parameterArray, q_reorient, origin, ntor, tlist, vt, &ligand, debug );
  1749.  
  1750.         (void) fflush(logFile);
  1751.         break;
  1752.  
  1753.  
  1754. //______________________________________________________________________________
  1755.  
  1756.     case DPF_TRAN0:
  1757.         /*
  1758.         **  tran0
  1759.         **  Initial_translation,
  1760.         */
  1761.         (void) sscanf( line, "%*s %s", param[0]);
  1762.         for (i=0; i<6; i++) {
  1763.             param[0][i] = (char)tolower( (int)param[0][i] );
  1764.         }
  1765.         if (equal(param[0],"random",6)) {
  1766.             B_RandomTran0 = TRUE;
  1767.             ligand.S.T.x = sInit.T.x = random_range( info->lo[X], info->hi[X] );
  1768.             ligand.S.T.y = sInit.T.y = random_range( info->lo[Y], info->hi[Y] );
  1769.             ligand.S.T.z = sInit.T.z = random_range( info->lo[Z], info->hi[Z] );
  1770.         } else {
  1771.             B_RandomTran0 = FALSE;
  1772.             (void) sscanf( line,"%*s %lf %lf %lf", &(sInit.T.x), &(sInit.T.y), &(sInit.T.z));
  1773.             ligand.S.T.x = sInit.T.x;
  1774.             ligand.S.T.y = sInit.T.y;
  1775.             ligand.S.T.z = sInit.T.z;
  1776.             if (outlev >= 0) {
  1777.                 pr( logFile, "Initial translation =\t\t\t(%.3f, %.3f, %.3f) Angstroms\n", sInit.T.x, sInit.T.y, sInit.T.z );
  1778.             }
  1779.             (void) fflush(logFile);
  1780.         }
  1781.         break;
  1782.  
  1783. //______________________________________________________________________________
  1784.  
  1785.     case DPF_AXISANGLE0:
  1786.     case DPF_QUATERNION0:
  1787.         /*
  1788.          * Handles both axisangle0 and quaternion0
  1789.          *
  1790.          *  axisangle0 1. 0. 0. 0.
  1791.          *  axisangle0 random
  1792.          *  ( quat0 <--- deprecated )
  1793.          *  Initial_quaternion, specified as an axis and angle
  1794.          *
  1795.          *  quaternion0 0. 0. 0. 1.
  1796.          *  quaternion0 random
  1797.          *  Initial_quaternion, specified as the four components (qx, qy, qz, qw)
  1798.          */
  1799.         {
  1800.         // Local Block...
  1801.         double a, b, c, d;
  1802.         (void) sscanf( line, "%*s %s", param[0]);
  1803.         for (i=0; i<6; i++) {
  1804.             param[0][i] = (char)tolower( (int)param[0][i] );
  1805.         }
  1806.         if (equal(param[0],"random",6)) {
  1807.             // Make a random initial quaternion,
  1808.             // and set the boolean B_RandomQuat0 to true,
  1809.             // so we can generate random quaternions in population-based methods.
  1810.             B_RandomQuat0 = TRUE;
  1811.             create_random_orientation( &(sInit.Q) );
  1812.             if (outlev >= 0) {
  1813.                 pr( logFile, "Each run will begin with a new, random initial orientation.\n");
  1814.             }
  1815.         } else {
  1816.             // Read in the user-defined axis-angle values for the initial quaternion
  1817.             // and set the boolean B_RandomQuat0 to false,
  1818.             B_RandomQuat0 = FALSE;
  1819.             (void) sscanf( line, "%*s %lf %lf %lf %lf", &a, &b, &c, &d);
  1820.             sInit.Q = (dpf_keyword == DPF_QUATERNION0) ?
  1821.                       quatComponentsToQuat(a,b,c,d) :
  1822.                       axisDegreeToQuat(a,b,c,d);
  1823.         }
  1824.         ligand.S.Q = sInit.Q;
  1825.         if (outlev >= 0) {
  1826.             if (dpf_keyword == DPF_QUATERNION0) {
  1827.                 pr( logFile, "Initial quaternion,  (x,y,z,w) =\t( %.3f, %.3f, %.3f, %.3f ),\n", sInit.Q.x, sInit.Q.y, sInit.Q.z, sInit.Q.w);
  1828.             } else {
  1829.                 pr( logFile, "Initial axis-angle,  (nx,ny,nz,ang) =\t( %.3f, %.3f, %.3f, %.1f deg ),\n", a, b, c, d );
  1830.                 pr( logFile, "Normalized axis,     (nx,ny,nz)     =\t( %.3f %.3f %.3f )\n", sInit.Q.nx, sInit.Q.ny, sInit.Q.nz );
  1831.             }
  1832. #ifdef DEBUG
  1833.             pr( logFile, "Initial Quaternion sInit.Q:\n\n");
  1834.             printQuat( logFile, sInit.Q );
  1835.             pr( logFile, "Initial Quaternion ligand.S.Q:\n\n");
  1836.             printQuat( logFile, ligand.S.Q );
  1837. #endif
  1838.         }
  1839.         (void) fflush(logFile);
  1840.         }
  1841.         break;
  1842.  
  1843. //______________________________________________________________________________
  1844.  
  1845.     case DPF_NDIHE:
  1846.         /*
  1847.         **  ndihe
  1848.         **  Number of dihedral angles to be specified by "dihe0"
  1849.         */
  1850.         (void) sscanf( line, "%*s %d", &ndihed );
  1851.         if ( nmol == 0 ) {
  1852.             if (outlev >= 0) {
  1853.                 pr( logFile, "Must specify a ligand PDBQT file, using the \"move\" command.\n");
  1854.             }
  1855.         } else {
  1856.             if (outlev >= 0) {
  1857.                 pr( logFile, "%s: WARNING!  The \"ndihe\" command is no longer supported.  The number of torsions in the PDBQT file(s) is the number that will be used (i.e. %d)\n", programname, ntor);
  1858.             }
  1859.             if ( ndihed != ntor ) {
  1860.                 pr( logFile, "%s: WARNING!  You requested %d torsions, but I found %d in PDBQT-file specifications.\n", programname, ndihed, ntor );
  1861.             } /* if */
  1862.         }
  1863.         (void) fflush(logFile);
  1864.         break;
  1865.  
  1866. //______________________________________________________________________________
  1867.  
  1868.     case DPF_DIHE0:
  1869.         /*
  1870.         **  dihe0
  1871.         **  Initial dihedral angles, input in degrees,
  1872.         */
  1873.         (void) sscanf( line, "%*s %s", param[0]);
  1874.         for (i=0; i<6; i++) {
  1875.             param[0][i] = (char)tolower( (int)param[0][i] );
  1876.         }
  1877.         if (equal(param[0],"random",6)) {
  1878.             B_RandomDihe0 = TRUE;
  1879.             sInit.ntor = nval = ntor;
  1880.             for ( i=0; i<nval; i++ ) {
  1881.                 sInit.tor[i] = random_range( -180.0, 180.0 );
  1882.             }
  1883.         } else {
  1884.             B_RandomDihe0 = FALSE;
  1885.             retval = (int)sscanf( line, torfmt, TOR_ARG_LIST );
  1886.             if (retval == 0) {
  1887.                 pr( logFile, "WARNING!  AutoDock could not read any torsions!\n" );
  1888.             } else if (retval == EOF) {
  1889.                 pr( logFile, "WARNING!  End of file encountered while reading dihe0 line\n");
  1890.             } else if (retval < ntor) {
  1891.                 pr( logFile, "WARNING!  Only %d initial torsion angles were detected on input line.\n",retval);
  1892.                 pr( logFile, "WARNING!  I am sorry, the number of torsions detected in the PDBQT files was %d torsions.\n", ntor);
  1893.             } else {
  1894.                 if (outlev >= 0) {
  1895.                     pr( logFile, "%d initial torsion angles were detected on input line.\n", retval );
  1896.                 }
  1897.             }
  1898.             nval = retval;
  1899.         }
  1900.         for ( i=0; i<nval; i++ ) {
  1901.             if (outlev >= 0) {
  1902.                 pr( logFile, "\tInitial torsion %2d = %7.2f deg\n", (i+1), sInit.tor[i] ); /* sInit.tor is in degrees */
  1903.                 /* Convert sInit.tor[i] into radians */
  1904.             }
  1905.             ligand.S.tor[i] = sInit.tor[i] = DegreesToRadians( sInit.tor[i] ); /* sInit.tor is now in radians  Added:05-01-95 */
  1906.         }
  1907.         (void) fflush(logFile);
  1908.         break;
  1909.  
  1910. //______________________________________________________________________________
  1911.  
  1912.     case DPF_TSTEP:
  1913.         /*
  1914.         **  tstep
  1915.         **  Translation_step,
  1916.         */
  1917.         retval = (int)sscanf( line, "%*s " FDFMT2, &trnStep0, &trnStepFinal );
  1918.         if (retval == 0) {
  1919.             pr( logFile, "WARNING!  AutoDock could not read any arguments!\n" );
  1920.         } else if (retval == EOF) {
  1921.             pr( logFile, "WARNING!  End of file encountered!\n");
  1922.         } else if (retval > 0) {
  1923.             pr( logFile, "Initial cycle, maximum translation step = +/- %-.1f Angstroms\n", trnStep0);
  1924.         }
  1925.         if (retval == 2) {
  1926.             B_CalcTrnRF = TRUE;
  1927.             if (outlev >= 0) {
  1928.                 pr( logFile, "Final cycle,   maximum translation step = +/- %-.1f Angstroms\n", trnStepFinal);
  1929.                 pr( logFile, "Reduction factor will be calculated when number of cycles has been read in.\n");
  1930.             }
  1931.         }
  1932.         (void) fflush(logFile);
  1933.         break;
  1934.  
  1935. //______________________________________________________________________________
  1936.  
  1937.     case DPF_QSTEP:
  1938.         /*
  1939.         **  qstep
  1940.         **  Quaternion_step,
  1941.         */
  1942.         retval = (int)sscanf( line, "%*s " FDFMT2, &qtwStep0, &qtwStepFinal );
  1943.         if (retval == 0) {
  1944.             pr( logFile, "WARNING!  AutoDock could not read any arguments!\n" );
  1945.         } else if (retval == EOF) {
  1946.             pr( logFile, "WARNING!  End of file encountered!\n");
  1947.         } else if (retval > 0) {
  1948.             if (outlev >= 0) {
  1949.                 pr( logFile, "Initial cycle, maximum quaternion angle step = +/- %-.1f deg\n", qtwStep0);
  1950.             }
  1951.             /* convert to radians */
  1952.             qtwStep0 = DegreesToRadians( qtwStep0 );
  1953.         }
  1954.         if (retval == 2) {
  1955.             B_CalcQtwRF = TRUE;
  1956.             if (outlev >= 0) {
  1957.                 pr( logFile, "Final cycle,   maximum quaternion angle step = +/- %-.1f deg\n", qtwStepFinal);
  1958.                 pr( logFile, "Reduction factor will be calculated when number of cycles has been read in.\n");
  1959.             }
  1960.             /* convert to radians */
  1961.             qtwStepFinal = DegreesToRadians( qtwStepFinal );
  1962.         }
  1963.         (void) fflush(logFile);
  1964.         break;
  1965.  
  1966. //______________________________________________________________________________
  1967.  
  1968.     case DPF_DSTEP:
  1969.         /*
  1970.         **  dstep
  1971.         **  Torsion_step,
  1972.         */
  1973.         retval = (int)sscanf( line, "%*s " FDFMT2, &torStep0, &torStepFinal );
  1974.         if (retval == 0) {
  1975.             pr( logFile, "WARNING!  AutoDock could not read any arguments!\n" );
  1976.         } else if (retval == EOF) {
  1977.             pr( logFile, "WARNING!  End of file encountered!\n");
  1978.         } else if (retval > 0) {
  1979.             if (outlev >= 0) {
  1980.                 pr( logFile, "Initial cycle, maximum torsion angle step = +/- %-.1f deg\n", torStep0);
  1981.             }
  1982.             /* convert to radians */
  1983.             torStep0 = DegreesToRadians( torStep0 );
  1984.         }
  1985.         if (retval == 2) {
  1986.             B_CalcTorRF = TRUE;
  1987.             if (outlev >= 0) {
  1988.                 pr( logFile, "Final cycle,   maximum torsion angle step = +/- %-.1f deg\n", torStepFinal);
  1989.                 pr( logFile, "Reduction factor will be calculated when number of cycles has been read in.\n");
  1990.             }
  1991.             /* convert to radians */
  1992.             torStepFinal = DegreesToRadians( torStepFinal );
  1993.         }
  1994.         (void) fflush(logFile);
  1995.         break;
  1996.  
  1997. //______________________________________________________________________________
  1998.  
  1999.     case DPF_TRNRF:
  2000.         /*
  2001.         **  trnrf
  2002.         **  Translation reduction factor,
  2003.         */
  2004.         (void) sscanf( line, "%*s " FDFMT, &trnFac );
  2005.         if (outlev >= 0) {
  2006.             pr( logFile, "Reduction factor for translations =\t%-.3f /cycle\n", trnFac );
  2007.         }
  2008.         B_trnReduc = (trnFac != 1.);
  2009.         (void) fflush(logFile);
  2010.         break;
  2011.  
  2012. //______________________________________________________________________________
  2013.  
  2014.     case DPF_QUARF:
  2015.         /*
  2016.         **  quarf
  2017.         **  Quaternion reduction factor,
  2018.         */
  2019.         (void) sscanf( line, "%*s " FDFMT, &qtwFac );
  2020.         if (outlev >= 0) {
  2021.             pr( logFile, "Reduction factor for quaternion angle =\t%-.3f /cycle\n", qtwFac );
  2022.         }
  2023.         B_qtwReduc = (qtwFac != 1.);
  2024.         (void) fflush(logFile);
  2025.         break;
  2026.  
  2027. //______________________________________________________________________________
  2028.  
  2029.     case DPF_DIHRF:
  2030.         /*
  2031.         **  dihrf
  2032.         **  Torsion reduction factor,
  2033.         */
  2034.         (void) sscanf( line, "%*s " FDFMT, &torFac );
  2035.         if (outlev >= 0) {
  2036.             pr( logFile, "Reduction factor for torsion angles =\t%-.3f /cycle\n", torFac );
  2037.         }
  2038.         B_torReduc = (torFac != 1.);
  2039.         (void) fflush(logFile);
  2040.         break;
  2041.  
  2042. //______________________________________________________________________________
  2043.  
  2044.     case DPF_FLEX:
  2045.         /*
  2046.         **  flex
  2047.         **  Flexible side-chains, cannot translate:
  2048.         */
  2049.         nmol++;
  2050.         nres++;
  2051.         (void) fflush(logFile);
  2052.         break;
  2053.  
  2054. //______________________________________________________________________________
  2055.  
  2056.     case DPF_INTNBP_REQM_EPS:
  2057.         /*
  2058.         **  intnbp_r_eps
  2059.         **  Read internal energy parameters:
  2060.         **  Lennard-Jones and Hydrogen Bond Potentials,
  2061.         **  Using epsilon and r-equilibrium values...
  2062.         */
  2063.         (void) sscanf( line, "%*s " FDFMT2 " %d %d", &Rij, &epsij, &xA, &xB );
  2064.         /* check that the Rij is reasonable */
  2065.         if ((Rij < RIJ_MIN) || (Rij > RIJ_MAX)) {
  2066.             (void) fprintf( logFile,
  2067.             "WARNING: pairwise distance, Rij, %.2f, is not a very reasonable value for the equilibrium separation of two atoms! (%.2f Angstroms <= Rij <= %.2f Angstroms)\n\n", Rij, RIJ_MIN, RIJ_MAX);
  2068.             (void) fprintf( logFile, "Perhaps you meant to use \"intnbp_coeffs\" instead of \"intnbp_r_eps\"?\n\n");
  2069.             /* GMM COMMENTED OUT FOR DAVE GOODSELL, MUTABLE ATOMS
  2070.              * exit(-1); */
  2071.         }
  2072.         /* check that the epsij is reasonable */
  2073.         if ((epsij < EPSIJ_MIN) || (epsij > EPSIJ_MAX)) {
  2074.             (void) fprintf( logFile,
  2075.             "WARNING: well-depth, epsilon_ij, %.2f, is not a very reasonable value for the equilibrium potential energy of two atoms! (%.2f kcal/mol <= epsilon_ij <= %.2f kcal/mol)\n\n", epsij, EPSIJ_MIN, EPSIJ_MAX);
  2076.             (void) fprintf( logFile, "Perhaps you meant to use \"intnbp_coeffs\" instead of \"intnbp_r_eps\"?\n\n");
  2077.             /* GMM COMMENTED OUT FOR DAVE GOODSELL, MUTABLE ATOMS
  2078.              * exit(-1); */
  2079.         }
  2080.         /* Defend against division by zero... */
  2081.         if (xA != xB) {
  2082.             // Calculate the coefficients from Rij and epsij
  2083.             cA = (tmpconst = epsij / (Real)(xA - xB)) * pow( (double)Rij, (double)xA ) * (Real)xB;
  2084.             cB = tmpconst * pow( (double)Rij, (double)xB ) * (Real)xA;
  2085.             pr(logFile, "\nCalculating internal non-bonded interaction energies for docking calculation;\n");
  2086.             intnbtable( &B_havenbp, a1, a2, info, cA, cB, xA, xB, AD4.coeff_desolv, sigma, ad_energy_tables, BOUND_CALCULATION );
  2087.             pr(logFile, "\nCalculating internal non-bonded interaction energies for unbound conformation calculation;\n");
  2088.             intnbtable( &B_havenbp, a1, a2, info, cA_unbound, cB_unbound, xA_unbound, xB_unbound, AD4.coeff_desolv, sigma, unbound_energy_tables, UNBOUND_CALCULATION );
  2089.             // Increment the atom type numbers, a1 and a2, for the internal non-bond table
  2090.             a2++;
  2091.             if (a2 >= info->num_atom_types) {
  2092.                 a1++;
  2093.                 a2 = a1;
  2094.             }
  2095.         } else {
  2096.             pr(logFile,"WARNING: Exponents must be different, to avoid division by zero!\n\tAborting...\n");
  2097.             exit(-1);
  2098.         }
  2099.         (void) fflush(logFile);
  2100.         break;
  2101.  
  2102. //______________________________________________________________________________
  2103.  
  2104.     case DPF_INTNBP_COEFFS:
  2105.         /*
  2106.         **  intnbp_coeffs
  2107.         **  Read internal energy parameters:
  2108.         **  Lennard-Jones and Hydrogen Bond Potentials,
  2109.         **  Using coefficients...
  2110.         */
  2111.         (void) sscanf( line, "%*s " FDFMT2 " %d %d", &cA, &cB, &xA, &xB );
  2112.  
  2113.         /* Defend against division by zero... */
  2114.         if (xA != xB) {
  2115.             pr(logFile, "\nCalculating internal non-bonded interaction energies for docking calculation;\n");
  2116.             intnbtable( &B_havenbp, a1, a2, info, cA, cB, xA, xB, AD4.coeff_desolv, sigma, ad_energy_tables, BOUND_CALCULATION );
  2117.             pr(logFile, "\nCalculating internal non-bonded interaction energies for unbound conformation calculation;\n");
  2118.             intnbtable( &B_havenbp, a1, a2, info, cA_unbound, cB_unbound, xA_unbound, xB_unbound, AD4.coeff_desolv, sigma, unbound_energy_tables, UNBOUND_CALCULATION );
  2119.             // Increment the atom type numbers, a1 and a2, for the internal non-bond table
  2120.             a2++;
  2121.             if (a2 >= info->num_atom_types) {
  2122.                 a1++;
  2123.                 a2 = a1;
  2124.             }
  2125.         } else {
  2126.             pr(logFile,"WARNING: Exponents must be different. Aborting...\n");
  2127.             exit(-1);
  2128.         }
  2129.  
  2130.         (void) fflush(logFile);
  2131.         break;
  2132.  
  2133. //______________________________________________________________________________
  2134.  
  2135.     case DPF_UNBOUND_INTNBP_COEFFS:
  2136.         /*
  2137.         **  unbound_intnbp_coeffs
  2138.         **  Read internal energy parameters for unbound extended state calculation:
  2139.         */
  2140.         (void) sscanf( line, "%*s " FDFMT2 " %d %d", &cA_unbound, &cB_unbound, &xA_unbound, &xB_unbound );
  2141.  
  2142.         pr(logFile, "\nSetting the internal non-bonded interaction energy parameters for the\nunbound docking calculation, E = %.1f / r^%d - %.1f / r^%d\n\n", cA_unbound, xA_unbound, cB_unbound, xB_unbound);
  2143.         (void) fflush(logFile);
  2144.         break;
  2145.  
  2146. //______________________________________________________________________________
  2147.  
  2148.     case DPF_RT0:
  2149.         /*
  2150.         **  rt0
  2151.         **  Initial Temperature,
  2152.         */
  2153.         (void) sscanf( line, "%*s " FDFMT, &RT0 );
  2154.         if (RT0 <= 0.) {
  2155.             pr( logFile, "\nWARNING!  Negative temperatures not allowed! Will default to RT = 616 cal mol.\n" );
  2156.             RT0 = 616.0;
  2157.         }
  2158.         if (outlev >= 0) {
  2159.             pr( logFile, "\n\t\tTEMPERATURE SCHEDULE INFORMATION\n" );
  2160.             pr( logFile, "\t\t________________________________\n\n" );
  2161.             pr( logFile, "               -1 -1                 -1 -1\n" );
  2162.             pr( logFile, "R = %5.3f J mol  K    = %5.3f cal mol  K  \n\n", RJ, Rcal );
  2163.             pr( logFile, "                                        -1\n" );
  2164.             pr( logFile, "Initial R*Temperature = %8.2f cal mol\n", RT0 );
  2165.             pr( logFile, "      (=> Temperature = %8.2f K\n", RT0/Rcal );
  2166.             pr( logFile, "                   or = %8.2f C)\n\n", RT0/Rcal - T0K );
  2167.         }
  2168.         break;
  2169.  
  2170. //______________________________________________________________________________
  2171.  
  2172.     case DPF_RTRF:
  2173.         /*
  2174.         **  rtrf
  2175.         **  Temperature reduction factor,
  2176.         */
  2177.         (void) sscanf( line, "%*s " FDFMT, &RTFac);
  2178.         if (outlev >= 0) {
  2179.             pr( logFile, "R*Temperature reduction factor = %8.2f\t/cycle\n", RTFac );
  2180.         }
  2181.         if (RTFac >= 1.) {
  2182.             stop("Cooling is impossible with a reduction\n\tfactor greater than or equal to 1.0!" );
  2183.             exit( -1 );
  2184.         } else if (RTFac == 0.0 ) {
  2185.             stop("Cooling is impossible with a ZERO reduction factor!" );
  2186.             exit( -1 );
  2187.         } else if (RTFac < 0.0 ) {
  2188.             stop("Cooling is impossible with a NEGATIVE reduction factor!" );
  2189.             exit( -1 );
  2190.         }
  2191.         (void) fflush(logFile);
  2192.         B_tempChange = ( RTFac != 1.0 );
  2193.         break;
  2194.  
  2195. //______________________________________________________________________________
  2196.  
  2197.     case DPF_RUNS:
  2198.         /*
  2199.         **  runs
  2200.         **  Number of docking runs,
  2201.         */
  2202.         (void) sscanf( line, "%*s %d", &nruns );
  2203.         if ( nruns > MAX_RUNS ) {
  2204.             prStr( error_message, "%s:  ERROR: %d runs were requested, but AutoDock is only dimensioned for %d.\nChange \"MAX_RUNS\" in \"constants.h\".", programname, nruns, MAX_RUNS);
  2205.             stop( error_message );
  2206.             exit( -1 );
  2207.         }
  2208.         pr( logFile, "Number of runs =\t\t\t\t%8d run%c\n", nruns, (nruns > 1)?'s':' ');
  2209.         (void) fflush(logFile);
  2210.         break;
  2211.  
  2212. //______________________________________________________________________________
  2213.  
  2214.     case DPF_CYCLES:
  2215.         /*
  2216.         **  cycles
  2217.         **  Number of constant temperature SA cycles,
  2218.         */
  2219.         (void) sscanf( line, "%*s %d", &ncycles );
  2220.         if (ncycles < 0) {
  2221.             pr( logFile, "WARNING!  Negative number of cycles found!  Using default value.\n");
  2222.             ncycles = 50;
  2223.         }
  2224.         pr( logFile, "Maximum number of cycles =\t\t\t%8d cycles\n\n", ncycles);
  2225.         if (B_linear_schedule) {
  2226.             RTreduc = RT0 / ncycles;
  2227.             if (outlev >= 0) {
  2228.                 pr( logFile, "\nA linear temperature reduction schedule was requested...\n" );
  2229.                 pr( logFile, "Annealing temperature will be reduced by %.3f cal mol per cycle.\n\n", RTreduc );
  2230.             }
  2231.         }
  2232.         (void) fflush(logFile);
  2233.         break;
  2234.  
  2235. //______________________________________________________________________________
  2236.  
  2237.     case DPF_ACCS:
  2238.         /*
  2239.         **  accs
  2240.         **  Maximum number of steps accepted,
  2241.         */
  2242.         (void) sscanf( line, "%*s %d", &naccmax );
  2243.         if (naccmax < 0) {
  2244.             naccmax = 100;
  2245.             pr( logFile, "WARNING!  Negative number of accepted moves found!  Using default value.\n");
  2246.         }
  2247.         if (outlev >= 0) {
  2248.             pr( logFile, "Maximum number accepted per cycle =\t\t%8d steps\n", naccmax);
  2249.         }
  2250.         if (nrejmax != 0) {
  2251.             nstepmax = naccmax + nrejmax;
  2252.             if (outlev >= 0) {
  2253.                 pr( logFile, "                                           \t_________\n" );
  2254.                 pr( logFile, "Maximum possible number of steps per cycle =\t%8d\tsteps\n\n", nstepmax);
  2255.             }
  2256.         }
  2257.         (void) fflush(logFile);
  2258.         break;
  2259.  
  2260. //______________________________________________________________________________
  2261.  
  2262.     case DPF_REJS:
  2263.         /*
  2264.         **  rejs
  2265.         **  Maximum number of steps rejected,
  2266.         */
  2267.         (void) sscanf( line, "%*s %d", &nrejmax );
  2268.         if (nrejmax < 0) {
  2269.             nrejmax = 100;
  2270.             pr( logFile, "WARNING!  Negative number of rejected moves found!  Using default value.\n");
  2271.         }
  2272.         if (outlev >= 0) {
  2273.             pr( logFile, "Maximum number rejected per cycle =\t\t%8d steps\n", nrejmax);
  2274.         }
  2275.         if (naccmax != 0) {
  2276.             nstepmax = naccmax + nrejmax;
  2277.             if (outlev >= 0) {
  2278.                 pr( logFile, "                                           \t_________\n" );
  2279.                 pr( logFile, "Maximum possible number of steps per cycle =\t%8d steps\n\n", nstepmax);
  2280.             }
  2281.         }
  2282.         (void) fflush(logFile);
  2283.         break;
  2284.  
  2285. //______________________________________________________________________________
  2286.  
  2287.     case DPF_SELECT:
  2288.         /*
  2289.         **  select
  2290.         **  Select either minimum or last state from previous cycle,
  2291.         */
  2292.         (void) sscanf( line, "%*s %1s", &selminpar );
  2293.         B_selectmin = (selminpar == 'm');
  2294.         if ( B_selectmin ) {
  2295.             if (outlev >= 0) {
  2296.                 pr( logFile, "%s will begin each new cycle\nwith the state of minimum energy from the previous annealing cycle.\n", programname);
  2297.             }
  2298.         } else {
  2299.             if (outlev >= 0) {
  2300.                 pr( logFile, "%s will begin each new cycle\nwith the last state from the previous annealing cycle.\n", programname);
  2301.             }
  2302.         }
  2303.         (void) fflush(logFile);
  2304.         break;
  2305.  
  2306. //______________________________________________________________________________
  2307.  
  2308.     case DPF_RMSTOL:
  2309.         /*
  2310.         **  rmstol
  2311.         **  Cluster tolerance,
  2312.         */
  2313.         (void) sscanf( line, "%*s " FDFMT, &clus_rms_tol);
  2314.         if (outlev >= 0) {
  2315.             pr( logFile, "Maximum RMS tolerance for conformational cluster analysis = %.2f Angstroms\n", clus_rms_tol);
  2316.         }
  2317.         (void) fflush(logFile);
  2318.         break;
  2319.  
  2320. //______________________________________________________________________________
  2321.  
  2322.     case DPF_RMSREF:
  2323.         /*
  2324.         **  rmsref
  2325.         **  RMS Reference Coordinates:
  2326.         */
  2327.         (void) sscanf( line, "%*s %s", FN_rms_ref_crds);
  2328.         if (outlev >= 0) {
  2329.             pr( logFile, "RMS reference coordinates will taken from \"%s\"\n", FN_rms_ref_crds );
  2330.         }
  2331.         (void) fflush(logFile);
  2332.         break;
  2333.  
  2334. //______________________________________________________________________________
  2335.  
  2336.     case DPF_RMSATOMS:
  2337.         /*
  2338.         **  rmsatoms ligand_only
  2339.         **  rmsatoms all
  2340.         **
  2341.         **  Set the atoms to compute the RMSD values for cluster analysis
  2342.         **  either "ligand_only" (the default) or "all" moving atoms (ligand + receptor)
  2343.         */
  2344.         retval = sscanf( line, "%*s %s", rms_atoms_cmd);
  2345.         if (retval != 1) {
  2346.             pr( logFile, "%s:  ERROR: please specify an argument (either \"ligand_only\" or \"all\").  By default, only the ligand atoms will be used for the cluster analysis.\n", programname );
  2347.             B_rms_atoms_ligand_only = TRUE;  // cluster on the ligand atoms only
  2348.         } else {
  2349.             if ( strncmp( rms_atoms_cmd, "ligand_only", 11 ) == 0 ) {
  2350.                 if (outlev >= 0) {
  2351.                     pr( logFile, "RMS clustering will be performed on the ligand atoms only.\n" );
  2352.                 }
  2353.                 B_rms_atoms_ligand_only = TRUE;  // cluster on the ligand atoms only
  2354.             } else if ( strncmp( rms_atoms_cmd, "all", 3 ) == 0 ) {
  2355.                 if (outlev >= 0) {
  2356.                     pr( logFile, "RMS clustering will be performed on the moving atoms of the receptor plus all the ligand atoms.\n" );
  2357.                 }
  2358.                 B_rms_atoms_ligand_only = FALSE;  // cluster on the ligand atoms plus moving receptor atoms
  2359.             } else {
  2360.                 if (outlev >= 0) {
  2361.                     pr( logFile, "RMS clustering will be performed on the ligand atoms only.\n" );
  2362.                 }
  2363.                 B_rms_atoms_ligand_only = TRUE;  // cluster on the ligand atoms only
  2364.             }
  2365.         }
  2366.         (void) fflush(logFile);
  2367.         break;
  2368.  
  2369. //______________________________________________________________________________
  2370.  
  2371.     case DPF_TRJFRQ:
  2372.         /*
  2373.         **  trjfrq
  2374.         **  Trajectory frequency,
  2375.         */
  2376.         (void) sscanf( line, "%*s %d", &trj_freq);
  2377.         B_write_trj = (trj_freq > 0);
  2378.         if (outlev >= 0) {
  2379.             pr( logFile, UnderLine );
  2380.             pr( logFile, "\t\tTRAJECTORY INFORMATION\n" );
  2381.             pr( logFile, "\t\t______________________\n\n\n" );
  2382.         }
  2383.         if (B_write_trj) {
  2384.             if (outlev >= 0) {
  2385.                 pr( logFile, "Output frequency for trajectory frames =\tevery %d step%s\n", trj_freq, (trj_freq > 1)?"s.":"." );
  2386.             }
  2387.         } else {
  2388.             if (outlev >= 0) {
  2389.                 pr( logFile, "No trajectory of states will be written.\n\n" );
  2390.                 pr( logFile, "Subsequent \"trjbeg\", \"trjend\", \"trjout\" and \"trjsel\" parameters will be ignored.\n\n" );
  2391.             }
  2392.         }
  2393.         (void) fflush(logFile);
  2394.         break;
  2395.  
  2396. //______________________________________________________________________________
  2397.  
  2398.     case DPF_TRJBEG:
  2399.         /*
  2400.         **  trjbeg
  2401.         **  Trajectory begin cycle,
  2402.         */
  2403.         (void) sscanf( line, "%*s %d", &trj_begin_cyc );
  2404.         if (outlev >= 0) {
  2405.             pr( logFile, "Begin outputting trajectory of states at cycle:\t%d\n", trj_begin_cyc );
  2406.         }
  2407.         if (trj_begin_cyc < 0) {
  2408.             trj_begin_cyc = 0;
  2409.         } else if (trj_begin_cyc > ncycles) {
  2410.             trj_begin_cyc = trj_end_cyc = ncycles;
  2411.         }
  2412.         --trj_begin_cyc;
  2413.         (void) fflush(logFile);
  2414.         break;
  2415.  
  2416. //______________________________________________________________________________
  2417.  
  2418.     case DPF_TRJEND:
  2419.         /*
  2420.         **  trjend
  2421.         **  Trajectory end cycle,
  2422.         */
  2423.         (void) sscanf( line, "%*s %d", &trj_end_cyc );
  2424.         if (outlev >= 0) {
  2425.             pr( logFile, "Cease outputting trajectory of states at cycle:\t%d\n", trj_end_cyc );
  2426.         }
  2427.         if (trj_end_cyc > ncycles) {
  2428.             trj_end_cyc = ncycles;
  2429.         } else if (trj_end_cyc < 0) {
  2430.             trj_end_cyc = 1;
  2431.         }
  2432.         --trj_end_cyc;
  2433.         (void) fflush(logFile);
  2434.         break;
  2435.  
  2436. //______________________________________________________________________________
  2437.  
  2438.     case DPF_TRJOUT:
  2439.         /*
  2440.         **  trjout
  2441.         **  Trajectory file,
  2442.         */
  2443.         (void) sscanf( line, "%*s %s", FN_trj );
  2444.         if (outlev >= 0) {
  2445.             pr( logFile, "\nWrite trajectory of state variables to file: \"%s\"\n", FN_trj);
  2446.         }
  2447.         (void) fflush(logFile);
  2448.         break;
  2449.  
  2450. //______________________________________________________________________________
  2451.  
  2452.     case DPF_TRJSEL:
  2453.         /*
  2454.         **  trjsel
  2455.         **  Trajectory select,
  2456.         */
  2457.         (void) sscanf( line, "%*s %1s", &out_acc_rej );
  2458.         B_acconly = (out_acc_rej == 'A');
  2459.         B_either  = (out_acc_rej == 'E');
  2460.         if (B_acconly) {
  2461.             if (outlev >= 0) {
  2462.                 pr( logFile, "Output *accepted* states only.\n" );
  2463.             }
  2464.         } else if (B_either) {
  2465.             if (outlev >= 0) {
  2466.                 pr( logFile, "Output *either* accepted or rejected states.\n" );
  2467.             }
  2468.         } else {
  2469.             pr( logFile, "WARNING: Missing or unknown accepted/rejected output flag.\n" );
  2470.         }
  2471.         (void) fflush(logFile);
  2472.         break;
  2473.  
  2474. //______________________________________________________________________________
  2475.  
  2476.     case DPF_EXTNRG:
  2477.         /*
  2478.         **  extnrg
  2479.         **  Wall Energy,
  2480.         */
  2481.         (void) sscanf( line, "%*s " FDFMT, &WallEnergy );
  2482.         if (outlev >= 0) {
  2483.             pr( logFile, "External grid energy (beyond grid map walls) = %.2f\n\n", WallEnergy );
  2484.         }
  2485.         (void) fflush(logFile);
  2486.         break;
  2487.  
  2488. //______________________________________________________________________________
  2489.  
  2490.     case DPF_CLUSTER:
  2491.         /*
  2492.         **  cluster
  2493.         **  Cluster mode,
  2494.         */
  2495.         (void) sscanf( line, "%*s %s", FN_clus );
  2496.         B_cluster_mode = TRUE;
  2497.         if (outlev >= 0) {
  2498.             pr( logFile, "Cluster mode is now set.\n\n" );
  2499.         }
  2500.         clmode( num_atom_types, clus_rms_tol,
  2501.                 hostnm, jobStart, tms_jobStart,
  2502.                 B_write_all_clusmem, FN_clus, crdpdb, lig_center,
  2503.                 B_symmetry_flag, FN_rms_ref_crds );
  2504.         (void) fflush(logFile);
  2505.         break;
  2506.  
  2507. //______________________________________________________________________________
  2508.  
  2509.     case DPF_CLUSALL:
  2510.         /*
  2511.         ** write_all_clusmem
  2512.         ** Write all cluster members...
  2513.         */
  2514.         B_write_all_clusmem = TRUE;
  2515.         if (outlev >= 0) {
  2516.             pr( logFile, "All members of each cluster will be written out after the clustering histogram.\n(This is instead of outputting just the lowest energy member in each.)\n\n" );
  2517.         }
  2518.         break;
  2519.  
  2520. //______________________________________________________________________________
  2521.  
  2522.     case DPF_RMSNOSYM:
  2523.         /*
  2524.         **  rmsnosym
  2525.         **  Calculate RMS values in the normal way,
  2526.         **  ignoring any atom-type equivalences...
  2527.         */
  2528.         B_symmetry_flag = FALSE;
  2529.         if (outlev >= 0) {
  2530.             pr( logFile, "Symmetry will be ignored in RMS calculations.\n\n" );
  2531.         }
  2532.         (void) fflush(logFile);
  2533.         break;
  2534.  
  2535. //______________________________________________________________________________
  2536.  
  2537.     case DPF_SCHEDLIN:
  2538.         /*
  2539.         **  linear_schedule
  2540.         **  Use a linear (arithmetic) temperature
  2541.         **  reduction schedule.  This is necessary for
  2542.         **  more accurate entropy estimations...
  2543.         */
  2544.         B_linear_schedule = TRUE;
  2545.         if (outlev >= 0) {
  2546.             pr( logFile, "A linear temperature reduction schedule will be used...\n\n" );
  2547.         }
  2548.         if (ncycles == -1) {
  2549.             pr( logFile, "\nWARNING!  Please specify the number of cycles first!\n\n" );
  2550.         } else {
  2551.             RTreduc = RT0 / ncycles;
  2552.             if (outlev >= 0) {
  2553.                 pr( logFile, "Annealing temperature will be reduced by %.3f cal mol per cycle.\n\n", RTreduc );
  2554.             }
  2555.         }
  2556.         (void) fflush(logFile);
  2557.         break;
  2558.  
  2559. //______________________________________________________________________________
  2560.  
  2561.     case DPF_WATCH:
  2562.         /*
  2563.         **  watch
  2564.         **  for watching a job's progress PDBQT file in AVS,
  2565.         */
  2566.         (void) sscanf( line, "%*s %s", FN_watch);
  2567.         if (B_write_trj) {
  2568.             pr(logFile,"\nAutoDock will create the watch-file \"%s\", for real-time monitoring of runs.\n\n", FN_watch);
  2569.             pr(logFile,"\nThe watch-file will be updated every %d moves, in accordance with the trajectory parameters..\n\n", trj_freq);
  2570.             B_watch = TRUE;
  2571.         } else {
  2572.             pr(logFile,"\nYou must set \"trjfrq\" to be greater than zero. No watch-file will be created.\n\n");
  2573.             B_watch = FALSE;
  2574.         }
  2575.         (void) fflush(logFile);
  2576.         break;
  2577.  
  2578. //______________________________________________________________________________
  2579.  
  2580.     case DPF_GAUSSTORCON:
  2581.     case DPF_HARDTORCON:
  2582.         /*
  2583.         ** "gausstorcon" Add Gaussian torsion contraints,
  2584.         ** "hardtorcon"  Add Hard torsion contraints,
  2585.         */
  2586.         (void) sscanf( line, "%*s %d " FDFMT2, &I_tor, &F_torPref, &F_torHWdth);
  2587.         if (I_tor <= 0) {
  2588.             pr( logFile, "\nTorsion IDs less than 1 (%d) are not allowed!\n\n", I_tor);
  2589.         } else if (I_tor > ntor) {
  2590.             pr( logFile, "\nRequested torsion ID (%d) is larger than the number of torsions found (%d)!\n\n", I_tor, ntor);
  2591.         } else { /* torsion-ID accepted */
  2592.             --I_tor;    /* Because humans start at 1, and C at 0... */
  2593.             if ( B_isTorConstrained[I_tor] == 0 ) {
  2594.  
  2595.                 if (dpf_keyword ==  DPF_GAUSSTORCON) {
  2596.                     B_isGaussTorCon = TRUE;
  2597.                     B_isTorConstrained[I_tor] = 1;
  2598.                     /*
  2599.                     ** Initialize... Torsion Energy Profile...
  2600.                     ** Set energies at every torsion division
  2601.                     ** to the user-defined (maximum) barrier energy,
  2602.                     */
  2603.                     for (US_tD = 0;  US_tD < NTORDIVS;  US_tD++) {
  2604.                         US_torProfile[I_tor][US_tD] = US_torBarrier;
  2605.                     }
  2606.                 } else {
  2607.                     /*
  2608.                     ** DPF_HARDTORCON
  2609.                     */
  2610.                     B_isTorConstrained[I_tor] = 2;
  2611.                 }
  2612.             }
  2613.             if (dpf_keyword ==  DPF_GAUSSTORCON) {
  2614.                 (void) strcpy( S_contype, " half-" );
  2615.             } else {
  2616.                 (void) strcpy( S_contype, " " );
  2617.             }
  2618.                 /*
  2619.             ** F_torPref ranges from -180.0 to +180.0 degrees...
  2620.             */
  2621.             F_torPref = WrpDeg(ModDeg(F_torPref));
  2622.             if (F_torHWdth < 0.) {
  2623.                 pr(logFile,"\nI'm sorry, negative%swidths (%.1f) are not allowed. I will use the default (%.1f) instead.\n\n", S_contype, F_torHWdth, DEFHWDTH);
  2624.                 F_torHWdth = DEFHWDTH;
  2625.             } else if (F_torHWdth > 90.) {
  2626.                 pr(logFile,"\nI'm sorry, your requested%swidth (%.1f) is too large. I will use the default (%.1f) instead.\n\n", S_contype, F_torHWdth, DEFHWDTH);
  2627.                 F_torHWdth = DEFHWDTH;
  2628.             }
  2629.             pr( logFile, "For torsion %d, Adding a constrained-torsion zone centered on %.1f degrees;\n%swidth = %.1f degrees.\n\n", 1+I_tor, F_torPref, S_contype, F_torHWdth);
  2630.  
  2631.             if (dpf_keyword == DPF_GAUSSTORCON) {
  2632.                 /*
  2633.                 ** Calculate the torsion energy profile;
  2634.                 ** combine this with previous profile without
  2635.                 ** losing any information.
  2636.                 */
  2637.                 for (F_A = F_A_from;  F_A <= F_A_to;  F_A += F_W) {
  2638.                     F_Aova = (F_A - F_torPref) / F_torHWdth;
  2639.                     US_energy = (unsigned short) (((Real)US_torBarrier) * (1.0 - exp(F_lnH * F_Aova*F_Aova)));
  2640.                     /*
  2641.                     ** if F_A(<-180.or>180), wrap to -180to180,
  2642.                     */
  2643.                     F_tor = WrpDeg(ModDeg(F_A));
  2644.                     /*
  2645.                     ** Convert from F_tor to US_tD
  2646.                     */
  2647.                     US_tD = (unsigned short) ((F_tor - F_hW + 180.)/F_W);
  2648.                     US_torProfile[I_tor][US_tD] = min(US_energy,US_torProfile[I_tor][US_tD]);
  2649.                 }/* for F_A */
  2650.                 /*
  2651.                 ** Ensure that the lowest point(s) in the profile are
  2652.                 ** zero...
  2653.                 */
  2654.                 US_min = TORBARMAX;
  2655.                 for (US_tD = 0;  US_tD < NTORDIVS;  US_tD++) {
  2656.                     US_min = min(US_min,US_torProfile[I_tor][US_tD]);
  2657.                 }
  2658.                 for (US_tD = 0;  US_tD < NTORDIVS;  US_tD++) {
  2659.                     US_torProfile[I_tor][US_tD] -= US_min;
  2660.                 }
  2661.             } else { /*DPF_HARDTORCON*/
  2662.  
  2663.                 iCon = N_con[I_tor] + 1;
  2664.                 if (iCon < MAX_TOR_CON) {
  2665.                     F_TorConRange[I_tor][N_con[I_tor]][LOWER] = F_torPref - 0.5* F_torHWdth;
  2666.                     F_TorConRange[I_tor][N_con[I_tor]][UPPER] = F_torPref + 0.5* F_torHWdth;
  2667.                     N_con[I_tor] = iCon;
  2668.                 } else {
  2669.                     pr(logFile,"\n\n I'm sorry, you can only have %d (=MAX_TOR_CON) torsion constraints.\nIf you need more, change the \"#define MAX_TOR_CON\" line in \"constants.h\".\n\n",MAX_TOR_CON);
  2670.                 }/* Still room to add another constraint. */
  2671.             } /*DPF_HARDTORCON*/
  2672.         }/* torsion-ID accepted */
  2673.         (void) fflush(logFile);
  2674.         break;
  2675.  
  2676. //______________________________________________________________________________
  2677.  
  2678.     case DPF_BARRIER:
  2679.         /*
  2680.         **  barrier
  2681.         **  Define torsion-barrier energy...
  2682.         */
  2683.         (void) sscanf( line, "%*s %d", &I_torBarrier);
  2684.         US_torBarrier = (unsigned short)I_torBarrier;
  2685.         US_torBarrier = min(US_torBarrier, TORBARMAX);
  2686.         pr(logFile,"\nTorsion barrier energy is set to %uhd\n\n", US_torBarrier);
  2687.         (void) fflush(logFile);
  2688.         break;
  2689.  
  2690. //______________________________________________________________________________
  2691.  
  2692.     case DPF_SHOWTORPEN:
  2693.         /*
  2694.         **  showtorpen
  2695.         **  Show torsion's penalty energy.
  2696.         */
  2697.         B_ShowTorE = TRUE;
  2698.         pr(logFile,"\nConstrained torsion penalty energies will be stored during docking, and output after each run\n\n");
  2699.         (void) fflush(logFile);
  2700.         break;
  2701.  
  2702. //______________________________________________________________________________
  2703.  
  2704.     case DPF_E0MAX:
  2705.         /*
  2706.         **  e0max
  2707.         **  Set maximum initial energy,
  2708.         */
  2709.         retval = sscanf( line, "%*s " FDFMT " %d", &e0max, &MaxRetries );
  2710.         if (retval == 0) {
  2711.             pr( logFile, "Could not read any arguments!\n" );
  2712.         } else if (retval == EOF) {
  2713.             pr( logFile, "End of file encountered!\n");
  2714.         } else if (retval == 1) {
  2715.             pr(logFile,"Using the default maximum number of retries for initialization, %d retries.\n\n", MaxRetries);
  2716.         } else if (retval == 2) {
  2717.             pr(logFile,"Using user-specified maximum number of retries for initialization, %d retries.\n\n", MaxRetries);
  2718.         }
  2719.         if (e0max < 0.) {
  2720.             e0max = 1000.0;
  2721.         }
  2722.         pr(logFile,"\nIf the initial energy is greater than e0max, %.3f,\nthen a new, random initial state will be created.\n\n",e0max);
  2723.         (void) fflush(logFile);
  2724.         break;
  2725.  
  2726. //______________________________________________________________________________
  2727.  
  2728.     case DPF_SIMANNEAL:
  2729.         /*
  2730.         ** simanneal
  2731.         */
  2732.         /*
  2733.         ** Calculate reduction factor based on initial and final step values,
  2734.         ** and number of cycles...
  2735.         */
  2736.         if (!command_mode) {
  2737.             ncycm1 = ncycles - 1;
  2738.             if (ncycm1 < 0) {
  2739.                 ncycm1 = 1;
  2740.             }
  2741.             if (B_CalcTrnRF) {
  2742.                 trnFac = RedFac(trnStep0, trnStepFinal, ncycm1);
  2743.                 pr( logFile, "Calculated reduction factor for translations     = %-.3f /cycle\n", trnFac);
  2744.                 B_trnReduc = (trnFac != 1.);
  2745.             }
  2746.             if (B_CalcQtwRF) {
  2747.                 qtwFac = RedFac(qtwStep0, qtwStepFinal, ncycm1);
  2748.                 pr( logFile, "Calculated reduction factor for quaternion angle = %-.3f /cycle\n", qtwFac );
  2749.                 B_qtwReduc = (qtwFac != 1.);
  2750.             }
  2751.             if (B_CalcTorRF) {
  2752.                 torFac    = RedFac(torStep0, torStepFinal, ncycm1);
  2753.                 pr( logFile, "Calculated reduction factor for torsion angles   = %-.3f /cycle\n", torFac );
  2754.                 B_torReduc = (torFac != 1.);
  2755.             }
  2756.             pr(logFile, "\n");
  2757.             /*
  2758.             ** Number of ligands read in...
  2759.             */
  2760.             if (nlig > 0) {
  2761.                 pr( logFile, "\nTotal number of ligands read in by the DPF \"move\" command = %d\n\n", nlig );
  2762.             }
  2763.             if (nres > 0) {
  2764.                 pr( logFile, "\nTotal number of residues read in by the DPF \"flex\" command = %d\n\n", nres );
  2765.             }
  2766.             pr(logFile, "\n");
  2767.             if (B_havenbp) {
  2768.                 nbe( info, ad_energy_tables, num_atom_types );
  2769.             }
  2770.             if (B_cluster_mode) {
  2771.                 clmode( num_atom_types, clus_rms_tol,
  2772.                         hostnm, jobStart, tms_jobStart,
  2773.                         B_write_all_clusmem, FN_clus, crdpdb, lig_center,
  2774.                         B_symmetry_flag, FN_rms_ref_crds );
  2775.             }
  2776.             for (j = 0; j < MAX_RUNS; j++) {
  2777.                 econf[j] = torsFreeEnergy - unbound_internal_FE;
  2778.             }
  2779.             /* ___________________________________________________________________
  2780.             **
  2781.             ** Begin the automated docking simulation,
  2782.             ** ___________________________________________________________________
  2783.             */
  2784.             simanneal( &nconf, Nnb, WallEnergy, atomstuff, charge, abs_charge, qsp_abs_charge, B_calcIntElec,
  2785.                         crd, crdpdb, dock_param_fn,
  2786.                         ad_energy_tables,
  2787.                         econf, B_either,
  2788.                         elec, emap,
  2789.                         ncycles, nruns, jobStart,
  2790.                         map,
  2791.                         naccmax, natom, nonbondlist, nrejmax, ntor1, ntor, outlev,
  2792.                         sInit, sHist,   qtwFac, B_qtwReduc, qtwStep0,
  2793.                         B_selectmin, FN_ligand, lig_center, RT0, B_tempChange, RTFac,
  2794.                         tms_jobStart, tlist, torFac, B_torReduc, torStep0,
  2795.                         FN_trj, trj_end_cyc, trj_begin_cyc, trj_freq, trnFac,
  2796.                         B_trnReduc, trnStep0, type, vt, B_write_trj,
  2797.                         B_constrain_dist, atomC1, atomC2, sqlower, squpper,
  2798.                         B_linear_schedule, RTreduc,
  2799.                         /*maxrad,*/
  2800.                         B_watch, FN_watch,
  2801.                         B_isGaussTorCon, US_torProfile, B_isTorConstrained,
  2802.                         B_ShowTorE, US_TorE, F_TorConRange, N_con,
  2803.                         B_RandomTran0, B_RandomQuat0, B_RandomDihe0,
  2804.                         e0max, torsFreeEnergy, MaxRetries, ligand_is_inhibitor,
  2805.                         ignore_inter,
  2806.                         B_include_1_4_interactions, scale_1_4,
  2807.                         parameterArray, unbound_internal_FE,
  2808.                         info, B_use_non_bond_cutoff,
  2809.                         B_have_flexible_residues,
  2810.                         PDBQT_record);
  2811.  
  2812.             (void) fflush(logFile);
  2813.         } else {
  2814.             (void)fprintf(logFile, "NOTE: Command mode has been set, so simulated annealing cannot be performed.\n\n");
  2815.             (void) fflush(logFile);
  2816.         }
  2817.         break;
  2818.  
  2819. //______________________________________________________________________________
  2820.  
  2821.     case DPF_SET_GA:
  2822.  
  2823.       if (GlobalSearchMethod != NULL) {
  2824.           pr(logFile, "Deleting the previous settings for the Genetic Algorithm.\n");
  2825.           (void) fflush(logFile);
  2826.           delete GlobalSearchMethod;
  2827.           GlobalSearchMethod = NULL;
  2828.       }
  2829.  
  2830.       if (debug > 0) {
  2831.           pr( logFile, "\n\tOutput every %u generations.\n", outputEveryNgens );
  2832.       }
  2833.       GlobalSearchMethod = new Genetic_Algorithm(e_mode, s_mode, c_mode, w_mode, elitism, c_rate, m_rate,
  2834.                                                  window_size, num_generations, outputEveryNgens );
  2835.       ((Genetic_Algorithm *)GlobalSearchMethod)->mutation_values( low, high, alpha, beta, trnStep0, qtwStep0, torStep0  );
  2836.       ((Genetic_Algorithm *)GlobalSearchMethod)->initialize(pop_size, 7+sInit.ntor);
  2837.  
  2838.       (void) fflush(logFile);
  2839.       break;
  2840. //______________________________________________________________________________
  2841.  
  2842.     case DPF_SET_SW1:
  2843.  
  2844.       if (LocalSearchMethod != NULL) {
  2845.           pr(logFile, "Deleting the previous settings for the local search Solis-Wets algorithm (SW1 object).\n");
  2846.           delete LocalSearchMethod;
  2847.           LocalSearchMethod = NULL;
  2848.       }
  2849.  
  2850.       pr(logFile, "Creating a new Local Search object using the Solis-Wets algorithm (SW1) with the current settings.\n\n");
  2851.       LocalSearchMethod = new Solis_Wets1(7+sInit.ntor, max_its, max_succ, max_fail, rho, lb_rho, 2.0, 0.5, search_freq);
  2852.  
  2853.       (void) fflush(logFile);
  2854.       break;
  2855. //______________________________________________________________________________
  2856.  
  2857.     case DPF_SET_PSW1:
  2858.  
  2859.       if (LocalSearchMethod != NULL) {
  2860.           pr(logFile, "Deleting the previous settings for the local search pseudo-Solis-Wets algorithm (pSW1 object).\n");
  2861.           delete LocalSearchMethod;
  2862.           LocalSearchMethod = NULL;
  2863.       }
  2864.  
  2865.       pr(logFile, "Creating a new Local Search object using the pseudo-Solis-Wets algorithm (pSW1) with the current settings.\n\n");
  2866.  
  2867.       //  Allocate space for the variable rho's
  2868.       rho_ptr = new Real[7+sInit.ntor];
  2869.       lb_rho_ptr = new Real[7+sInit.ntor];
  2870.  
  2871.       //  Initialize the rho's corresponding to the translation
  2872.       //  0,1,2   x,y,z
  2873.       //  3,4,5,6 qx,qy,qz,qw
  2874.       //  7,...   tor1
  2875.       for (j=0; j<3; j++) {
  2876.          // j=0,1,2
  2877.          rho_ptr[j] = trnStep0;
  2878.          // lb_rho_ptr[j] = trnStepFinal;
  2879. #define LB_RHO_FRACTION 0.10 // fraction of step size for Solis & Wets local search
  2880.          lb_rho_ptr[j] = LB_RHO_FRACTION * trnStep0;
  2881.       }
  2882.  
  2883.       //  Initialize the rho's corresponding to the quaterion
  2884.       for (; j<7; j++) {
  2885.          // j=3,4,5,6
  2886.          // rho_ptr[j] = qtwStep0;
  2887.          rho_ptr[j] = 1.0; // 4D normalized quaternion's components
  2888.          // lb_rho_ptr[j] = qtwStepFinal;
  2889.          // lb_rho_ptr[j] = LB_RHO_FRACTION * qtwStep0;
  2890.          lb_rho_ptr[j] = LB_RHO_FRACTION * 1.0;
  2891.       }
  2892.  
  2893.       //  Initialize the rho's corresponding to the torsions
  2894.       for (; j<7+sInit.ntor; j++) {
  2895.          // j=7,...
  2896.          rho_ptr[j] = torStep0;
  2897.          // lb_rho_ptr[j] = torStepFinal;
  2898.          lb_rho_ptr[j] = LB_RHO_FRACTION * torStep0;
  2899.       }
  2900.  
  2901.       LocalSearchMethod = new Pseudo_Solis_Wets1(7+sInit.ntor, max_its, max_succ, max_fail, 2.0, 0.5, search_freq, rho_ptr, lb_rho_ptr);
  2902.  
  2903.       (void) fflush(logFile);
  2904.       break;
  2905.  
  2906. //______________________________________________________________________________
  2907.  
  2908.     case DPF_SET_PATTERN:
  2909.  
  2910.       if (LocalSearchMethod != NULL) {
  2911.           pr(logFile, "Deleting the previous settings for the local search Pattern Search algorithm (PS object).\n");
  2912.           delete LocalSearchMethod;
  2913.           LocalSearchMethod = NULL;
  2914.       }
  2915.  
  2916.  
  2917.       pr(logFile, "Creating a new Local Search object using the Pattern Search algorithm (PS) with the current settings.\n\n");
  2918.       LocalSearchMethod = new Pattern_Search(7+sInit.ntor, max_succ, rho, lb_rho, 2.0, 0.5, search_freq);
  2919.  
  2920.       (void) fflush(logFile);
  2921.       break;
  2922.  
  2923. //______________________________________________________________________________
  2924.  
  2925.     case DPF_GALS:
  2926.         (void) fflush( logFile );
  2927.         /*
  2928.         ** Genetic Algorithm-Local search,  a.k.a. Lamarckian Genetic Algorithm
  2929.         */
  2930.         if (!command_mode) {
  2931.             (void) sscanf( line, "%*s %d", &nruns );
  2932.             if ( nruns > MAX_RUNS ) {
  2933.                 prStr( error_message, "%s:  ERROR: %d runs requested, but only dimensioned for %d.\nChange \"MAX_RUNS\" in \"constants.h\".", programname, nruns, MAX_RUNS);
  2934.                 stop( error_message );
  2935.                 exit( -1 );
  2936.             } else if ((GlobalSearchMethod==NULL)||(LocalSearchMethod==NULL)) {
  2937.                 prStr(error_message, "%s:  ERROR:  You must use \"set_ga\" to allocate both Global Optimization object AND Local Optimization object.\n", programname);
  2938.                 stop(error_message);
  2939.                 exit(-1);
  2940.             } else if (!B_found_elecmap) {
  2941.                 prStr(error_message, "%s:  ERROR:  no \"elecmap\" command has been specified!\n", programname);
  2942.                 stop(error_message);
  2943.                 exit(-1);
  2944.             } else if (!B_found_desolvmap) {
  2945.                 prStr(error_message, "%s:  ERROR:  no \"desolvmap\" command has been specified!\n", programname);
  2946.                 stop(error_message);
  2947.                 exit(-1);
  2948.             }
  2949.  
  2950.             pr( logFile, "Number of requested LGA dockings = %d run%c\n", nruns, (nruns > 1)?'s':' ');
  2951.  
  2952. #ifdef DEBUG
  2953.             pr( logFile, "\nAbout to call evaluate.setup(), sInit:\n\n");
  2954.             printState( logFile, sInit, 2 );
  2955. #endif
  2956.  
  2957.             evaluate.setup( crd, charge, abs_charge, qsp_abs_charge, type, natom, map,
  2958.                             elec, emap, nonbondlist, ad_energy_tables, Nnb,
  2959.                             B_calcIntElec, B_isGaussTorCon, B_isTorConstrained,
  2960.                             B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, crdreo, sInit, ligand,
  2961.                             ignore_inter,
  2962.                             B_include_1_4_interactions, scale_1_4,
  2963.                             parameterArray, unbound_internal_FE, info, B_use_non_bond_cutoff, B_have_flexible_residues);
  2964.  
  2965.             evaluate.compute_intermol_energy(TRUE);
  2966.  
  2967.             if(write_stateFile){
  2968.               fprintf(stateFile,"\t<run_requested>%d</run_requested>\n",nruns);
  2969.               fprintf(stateFile,"\t<runs>\n");
  2970.             }
  2971.             for (j=0; j<nruns; j++) {
  2972.                 j1=j+1;
  2973.  
  2974.                 (void) fprintf( logFile, "\n\n\tBEGINNING LAMARCKIAN GENETIC ALGORITHM DOCKING\n");
  2975.                 (void) fflush( logFile );
  2976.                 pr( logFile, "\nRun:\t%d / %d\n", j1, nruns );
  2977.  
  2978.                 pr(logFile, "Date:\t");
  2979.                 printdate( logFile, 2 );
  2980.                 (void) fflush( logFile );
  2981.  
  2982.                 gaStart = times( &tms_gaStart );
  2983.  
  2984.                 if (B_reorient_random == TRUE) {
  2985.                     // create a new random orientation before docking
  2986.                     create_random_orientation( &q_reorient );
  2987.                     // reorient the ligand
  2988.                     reorient( logFile, true_ligand_atoms, atomstuff, crdpdb, charge, type,
  2989.                               parameterArray, q_reorient, origin, ntor, tlist, vt, &ligand, debug );
  2990.                     // update the evaluate object
  2991.                     evaluate.update_crds( crdpdb, vt );
  2992.                 }
  2993.  
  2994.                 //  Can get rid of the following line
  2995.                 ((Genetic_Algorithm *)GlobalSearchMethod)->initialize(pop_size, 7+sInit.ntor);
  2996.  
  2997.                 // Reiterate output level...
  2998.                 pr(logFile, "Output level is set to %d.\n\n", outlev);
  2999.  
  3000.                 // Start Lamarckian GA run -- Bound simulation
  3001.                 sHist[nconf] = call_glss( GlobalSearchMethod, LocalSearchMethod,
  3002.                                           sInit,
  3003.                                           num_evals, pop_size,
  3004.                                           outlev,
  3005.                                           outputEveryNgens, &ligand,
  3006.                                           B_RandomTran0, B_RandomQuat0, B_RandomDihe0,
  3007.                                           info, FN_pop_file );
  3008.                 // State of best individual at end of GA-LS run is returned.
  3009.                 // Finished Lamarckian GA run
  3010.  
  3011.                 gaEnd = times( &tms_gaEnd );
  3012.                 pr( logFile, "\nRun completed;  time taken for this run:\n");
  3013.                 timesyshms( gaEnd - gaStart, &tms_gaStart, &tms_gaEnd );
  3014.                 pr( logFile, "\n");
  3015.                 printdate( logFile, 1 );
  3016.                 (void) fflush( logFile );
  3017.  
  3018.                 pr(logFile, "Total number of Energy Evaluations: %lu\n", evaluate.evals() );
  3019.                 pr(logFile, "Total number of Generations:        %u\n", ((Genetic_Algorithm *)GlobalSearchMethod)->num_generations());
  3020.  
  3021.                 pr( logFile, "\n\n\tFINAL LAMARCKIAN GENETIC ALGORITHM DOCKED STATE\n" );
  3022.                 pr( logFile,     "\t_______________________________________________\n\n\n" );
  3023.  
  3024.                 writePDBQT( j, seed,  FN_ligand, dock_param_fn, lig_center,
  3025.                             sHist[nconf], ntor, &eintra, &einter, natom, atomstuff,
  3026.                             crd, emap, elec,
  3027.                             charge, abs_charge, qsp_abs_charge,
  3028.                             ligand_is_inhibitor,
  3029.                             torsFreeEnergy,
  3030.                             vt, tlist, crdpdb, nonbondlist,
  3031.                             ad_energy_tables,
  3032.                             type, Nnb, B_calcIntElec,
  3033.                             map,
  3034.                             outlev, ignore_inter,
  3035.                             B_include_1_4_interactions, scale_1_4, parameterArray, unbound_internal_FE,
  3036.                             info, DOCKED, PDBQT_record, B_use_non_bond_cutoff, B_have_flexible_residues);
  3037.  
  3038.                 econf[nconf] = eintra + einter + torsFreeEnergy - unbound_internal_FE;
  3039.  
  3040.                 ++nconf;
  3041.  
  3042.                 pr( logFile, UnderLine );
  3043.             } // Next LGA run
  3044.             if(write_stateFile){
  3045.                fprintf(stateFile,"\t</runs>\n");
  3046.                (void) fflush(stateFile);
  3047.             }
  3048.             (void) fflush(logFile);
  3049.         } else {
  3050.             (void)fprintf(logFile, "NOTE: Command mode has been set.  Sorry, genetic algorithm-local search cannot be performed.\n\n");
  3051.         }
  3052.         break;
  3053.  
  3054. //______________________________________________________________________________
  3055.  
  3056.     case DPF_LS:
  3057.        // do_local_only
  3058.        (void) sscanf(line, "%*s %d", &nruns);
  3059.  
  3060.         if (!command_mode) {
  3061.             if ( nruns > MAX_RUNS ) {
  3062.  
  3063.                prStr( error_message, "%s:  ERROR: %d runs requested, but only dimensioned for %d.\nChange \"MAX_RUNS\" in \"constants.h\".", programname, nruns, MAX_RUNS);
  3064.                stop( error_message );
  3065.                exit( -1 );
  3066.  
  3067.            } else if (LocalSearchMethod==NULL) {
  3068.  
  3069.                prStr(error_message, "%s:  ERROR:  You must use \"set_sw1\", \"set_psw1\" or \"set_pattern\" to create a Local Optimization object.\n", programname);
  3070.                stop(error_message);
  3071.                exit(-1);
  3072.            } else if (!B_found_elecmap) {
  3073.                prStr(error_message, "%s:  ERROR:  no \"elecmap\" command has been specified!\n", programname);
  3074.                stop(error_message);
  3075.                exit(-1);
  3076.            } else if (!B_found_desolvmap) {
  3077.                prStr(error_message, "%s:  ERROR:  no \"desolvmap\" command has been specified!\n", programname);
  3078.                stop(error_message);
  3079.                exit(-1);
  3080.            }
  3081.  
  3082.  
  3083.            pr( logFile, "Number of Local Search (LS) only dockings = %d run%c\n", nruns, (nruns > 1)?'s':' ');
  3084.            evaluate.setup( crd, charge, abs_charge, qsp_abs_charge, type, natom, map,
  3085.                            elec, emap,
  3086.                            nonbondlist,
  3087.                            ad_energy_tables,
  3088.                            Nnb, B_calcIntElec, B_isGaussTorCon,B_isTorConstrained,
  3089.                            B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, crdreo, sInit, ligand,
  3090.                            ignore_inter,
  3091.                            B_include_1_4_interactions, scale_1_4,
  3092.                            parameterArray, unbound_internal_FE, info, B_use_non_bond_cutoff, B_have_flexible_residues);
  3093.  
  3094.             evaluate.compute_intermol_energy(TRUE);
  3095.  
  3096.            if(write_stateFile){
  3097.              fprintf(stateFile,"\t<run_requested>%d</run_requested>\n",nruns);
  3098.              fprintf(stateFile,"\t<runs>\n");
  3099.            }
  3100.  
  3101.            for (j=0; j<nruns; j++) {
  3102.  
  3103.                (void) fprintf( logFile, "\n\n\tBEGINNING SOLIS & WETS LOCAL SEARCH DOCKING\n");
  3104.                pr( logFile, "\nRun:\t%d / %d\n", j+1, nruns );
  3105.                (void) fflush( logFile );
  3106.  
  3107.                pr(logFile, "Date:\t");
  3108.                printdate( logFile, 2 );
  3109.                (void) fflush( logFile );
  3110.  
  3111.                gaStart = times( &tms_gaStart );
  3112.  
  3113.                sHist[nconf] = call_ls(LocalSearchMethod, sInit, pop_size, &ligand);
  3114.  
  3115.                pr(logFile, "There were %lu Energy Evaluations.\n", evaluate.evals());
  3116.  
  3117.                gaEnd = times( &tms_gaEnd );
  3118.                pr( logFile, "Time taken for this Local Search (LS) run:\n");
  3119.                timesyshms( gaEnd - gaStart, &tms_gaStart, &tms_gaEnd );
  3120.                pr( logFile, "\n");
  3121.                (void) fflush( logFile );
  3122.  
  3123.                pr( logFile, "\n\n\tFINAL LOCAL SEARCH DOCKED STATE\n" );
  3124.                pr( logFile,     "\t_______________________________\n\n\n" );
  3125.  
  3126.                writePDBQT( j, seed, FN_ligand, dock_param_fn, lig_center,
  3127.                            sHist[nconf], ntor, &eintra, &einter, natom, atomstuff,
  3128.                            crd, emap, elec,
  3129.                            charge, abs_charge, qsp_abs_charge,
  3130.                            ligand_is_inhibitor,
  3131.                            torsFreeEnergy,
  3132.                            vt, tlist, crdpdb, nonbondlist,
  3133.                            ad_energy_tables,
  3134.                            type, Nnb, B_calcIntElec,
  3135.                            map,
  3136.                            outlev,
  3137.                            ignore_inter,
  3138.                            B_include_1_4_interactions, scale_1_4, parameterArray, unbound_internal_FE,
  3139.                            info, DOCKED, PDBQT_record, B_use_non_bond_cutoff, B_have_flexible_residues);
  3140.  
  3141.                econf[nconf] = eintra + einter + torsFreeEnergy - unbound_internal_FE;
  3142.  
  3143.                ++nconf;
  3144.  
  3145.                pr( logFile, UnderLine );
  3146.  
  3147.            } // Next run
  3148.            if(write_stateFile){
  3149.              fprintf(stateFile,"\t</runs>\n");
  3150.              (void) fflush(stateFile);
  3151.            }
  3152.            (void) fflush(logFile);
  3153.        } else {
  3154.             (void)fprintf(logFile, "NOTE: Command mode has been set, so local search cannot be performed.\n\n");
  3155.        }
  3156.        break;
  3157.  
  3158. //______________________________________________________________________________
  3159.  
  3160.     case DPF_GS:
  3161.       (void) sscanf(line, "%*s %d", &nruns);
  3162.  
  3163.       if (!command_mode) {
  3164.           if (nruns>MAX_RUNS) {
  3165.  
  3166.               prStr(error_message, "%s:  ERROR:  %d runs requested, but only dimensioned for %d.\nChange \"MAX_RUNS\" in \"constants.h\".", programname, nruns, MAX_RUNS);
  3167.               stop(error_message);
  3168.               exit(-1);
  3169.  
  3170.           } else if (GlobalSearchMethod==NULL) {
  3171.               prStr(error_message, "%s:  ERROR:  You must use \"set_ga\" to allocate a Global Optimization object.\n", programname);
  3172.               stop(error_message);
  3173.               exit(-1);
  3174.           } else if (!B_found_elecmap) {
  3175.                prStr(error_message, "%s:  ERROR:  no \"elecmap\" command has been specified!\n", programname);
  3176.                stop(error_message);
  3177.                exit(-1);
  3178.            }
  3179.            else if (!B_found_desolvmap) {
  3180.                prStr(error_message, "%s:  ERROR:  no \"desolvmap\" command has been specified!\n", programname);
  3181.                stop(error_message);
  3182.                exit(-1);
  3183.            }
  3184.  
  3185.  
  3186.           pr(logFile, "Number of Genetic Algorithm (GA) only dockings = %d run%c\n", nruns, (nruns>1)?'s':' ');
  3187.  
  3188.  
  3189.           evaluate.setup( crd, charge, abs_charge, qsp_abs_charge, type, natom, map,
  3190.                           elec, emap,
  3191.                           nonbondlist,
  3192.                           ad_energy_tables,
  3193.                           Nnb, B_calcIntElec, B_isGaussTorCon,B_isTorConstrained,
  3194.                           B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, crdreo, sInit, ligand,
  3195.                           ignore_inter,
  3196.                           B_include_1_4_interactions, scale_1_4,
  3197.                           parameterArray, unbound_internal_FE, info, B_use_non_bond_cutoff, B_have_flexible_residues);
  3198.  
  3199.             evaluate.compute_intermol_energy(TRUE);
  3200.  
  3201.           if(write_stateFile){
  3202.             fprintf(stateFile,"\t<run_requested>%d</run_requested>\n",nruns);
  3203.             fprintf(stateFile,"\t<runs>\n");
  3204.           }
  3205.           for (j=0; j<nruns; j++) {
  3206.  
  3207.               fprintf( logFile, "\n\n\tBEGINNING GENETIC ALGORITHM DOCKING\n");
  3208.               pr(logFile, "\nDoing Genetic Algorithm run:  %d/%d.\n", j+1, nruns);
  3209.               (void) fflush( logFile );
  3210.  
  3211.               pr(logFile, "Date:\t");
  3212.               printdate( logFile, 2 );
  3213.               (void) fflush( logFile );
  3214.  
  3215.               gaStart = times(&tms_gaStart);
  3216.  
  3217.               sHist[nconf] = call_gs( GlobalSearchMethod, sInit, num_evals, pop_size,
  3218.                                       &ligand, outputEveryNgens, info);
  3219.  
  3220.               pr(logFile, "\nFinal docked state:\n");
  3221.               printState(logFile, sHist[nconf], 2);
  3222.  
  3223.               gaEnd = times(&tms_gaEnd);
  3224.               pr(logFile, "Time taken for this Genetic Algorithm (GA) run:\n");
  3225.               timesyshms(gaEnd-gaStart, &tms_gaStart, &tms_gaEnd);
  3226.               pr(logFile, "\n");
  3227.               (void) fflush(logFile);
  3228.  
  3229.               pr(logFile, "Total number of Energy Evaluations: %lu\n", evaluate.evals() );
  3230.               pr(logFile, "Total number of Generations:        %u\n", ((Genetic_Algorithm *)GlobalSearchMethod)->num_generations());
  3231.  
  3232.  
  3233.               pr( logFile, "\n\n\tFINAL GENETIC ALGORITHM DOCKED STATE\n" );
  3234.               pr( logFile,     "\t____________________________________\n\n\n" );
  3235.  
  3236.               writePDBQT( j, seed, FN_ligand, dock_param_fn, lig_center,
  3237.                           sHist[nconf], ntor, &eintra, &einter, natom, atomstuff,
  3238.                           crd, emap, elec,
  3239.                           charge, abs_charge, qsp_abs_charge,
  3240.                           ligand_is_inhibitor,
  3241.                           torsFreeEnergy,
  3242.                           vt, tlist, crdpdb, nonbondlist,
  3243.                           ad_energy_tables,
  3244.                           type, Nnb, B_calcIntElec,
  3245.                           map,
  3246.                           outlev,
  3247.                           ignore_inter,
  3248.                           B_include_1_4_interactions, scale_1_4, parameterArray, unbound_internal_FE,
  3249.                           info, DOCKED, PDBQT_record, B_use_non_bond_cutoff, B_have_flexible_residues);
  3250.  
  3251.               econf[nconf] = eintra + einter + torsFreeEnergy - unbound_internal_FE;
  3252.  
  3253.               ++nconf;
  3254.  
  3255.               pr( logFile, UnderLine );
  3256.  
  3257.           } // Next run
  3258.           if(write_stateFile){
  3259.             fprintf(stateFile,"\t</runs>\n");
  3260.             (void) fflush(stateFile);
  3261.           }
  3262.           (void) fflush(logFile);
  3263.       } else {
  3264.             (void)fprintf(logFile, "NOTE: Command mode has been set, so global search cannot be performed.\n\n");
  3265.       }
  3266.       break;
  3267.  
  3268. //______________________________________________________________________________
  3269.  
  3270.     case GA_pop_size:
  3271.        (void) sscanf(line, "%*s %u", &pop_size);
  3272.        pr(logFile, "A population of %u individuals will be used\n", pop_size);
  3273.         (void) fflush(logFile);
  3274.        break;
  3275.  
  3276. //______________________________________________________________________________
  3277.  
  3278.     case GA_num_generations:
  3279.        (void) sscanf(line, "%*s %u", &num_generations);
  3280.        pr(logFile, "The GA will run for at most %u generations.\n", num_generations);
  3281.         (void) fflush(logFile);
  3282.        break;
  3283.  
  3284. //______________________________________________________________________________
  3285.  
  3286.     case GA_num_evals:
  3287.        (void) sscanf(line, "%*s %u", &num_evals);
  3288.        pr(logFile, "There will be at most %u function evaluations used.\n", num_evals);
  3289.         (void) fflush(logFile);
  3290.        break;
  3291.  
  3292. //______________________________________________________________________________
  3293.  
  3294.     case GA_window_size:
  3295.        (void) sscanf(line, "%*s %d", &window_size);
  3296.        pr(logFile, "The GA's selection window is %d generations.\n", window_size);
  3297.         (void) fflush(logFile);
  3298.        break;
  3299.  
  3300. //______________________________________________________________________________
  3301.  
  3302.     case GA_low:
  3303.        (void) sscanf(line, "%*s %d", &low);
  3304.        pr(logFile, "Setting low to %d.\n", low);
  3305.         (void) fflush(logFile);
  3306.        break;
  3307.  
  3308. //______________________________________________________________________________
  3309.  
  3310.     case GA_high:
  3311.        (void) sscanf(line, "%*s %d", &high);
  3312.        pr(logFile, "Setting high to %d.\n", high);
  3313.         (void) fflush(logFile);
  3314.        break;
  3315.  
  3316. //______________________________________________________________________________
  3317.  
  3318.     case GA_elitism:
  3319.        (void) sscanf(line, "%*s %d", &elitism);
  3320.        pr(logFile, "The %d best will be preserved each generation.\n", elitism);
  3321.         (void) fflush(logFile);
  3322.        break;
  3323.  
  3324. //______________________________________________________________________________
  3325.  
  3326.     case GA_mutation_rate:
  3327.        (void) sscanf(line, "%*s " FDFMT, &m_rate);
  3328.        pr(logFile, "The mutation rate is %f.\n", m_rate);
  3329.         (void) fflush(logFile);
  3330.        break;
  3331.  
  3332. //______________________________________________________________________________
  3333.  
  3334.     case GA_crossover_rate:
  3335.        (void) sscanf(line, "%*s " FDFMT, &c_rate);
  3336.        pr(logFile, "The crossover rate is %f.\n", c_rate);
  3337.         (void) fflush(logFile);
  3338.        break;
  3339.  
  3340. //______________________________________________________________________________
  3341.  
  3342.     case GA_Cauchy_alpha:
  3343.        (void) sscanf(line, "%*s " FDFMT, &alpha);
  3344.        pr(logFile, "The alpha parameter (for the Cauchy distribution) is being set to %f.\n",
  3345.           alpha);
  3346.         (void) fflush(logFile);
  3347.        break;
  3348.  
  3349. //______________________________________________________________________________
  3350.  
  3351.     case GA_Cauchy_beta:
  3352.        (void) sscanf(line, "%*s " FDFMT, &beta);
  3353.        pr(logFile, "The beta parameter (for the Cauchy distribution) is being set to %f.\n",
  3354.           beta);
  3355.         (void) fflush(logFile);
  3356.        break;
  3357.  
  3358. //______________________________________________________________________________
  3359.  
  3360.     case SW_max_its:
  3361.        (void) sscanf(line, "%*s %u", &max_its);
  3362.        pr(logFile, "Solis & Wets algorithms will perform at most %u iterations.\n", max_its);
  3363.         (void) fflush(logFile);
  3364.        break;
  3365.  
  3366. //______________________________________________________________________________
  3367.  
  3368.     case SW_max_succ:
  3369.        (void) sscanf(line, "%*s %u", &max_succ);
  3370.        pr(logFile, "Solis & Wets algorithms expand rho every %u in a row successes.\n", max_succ);
  3371.         (void) fflush(logFile);
  3372.       break;
  3373.  
  3374. //______________________________________________________________________________
  3375.  
  3376.     case SW_max_fail:
  3377.        (void) sscanf(line, "%*s %u", &max_fail);
  3378.        pr(logFile, "Solis & Wets algorithms contract rho every %u in a row failures.\n", max_fail);
  3379.         (void) fflush(logFile);
  3380.       break;
  3381.  
  3382. //______________________________________________________________________________
  3383.  
  3384.     case SW_rho:
  3385.        (void) sscanf(line, "%*s " FDFMT, &rho);
  3386.        pr(logFile, "rho is set to %f.\n", rho);
  3387.         (void) fflush(logFile);
  3388.       break;
  3389.  
  3390. //______________________________________________________________________________
  3391.  
  3392.     case SW_lb_rho:
  3393.         (void) sscanf(line, "%*s " FDFMT, &lb_rho);
  3394.         pr(logFile, "rho will never get smaller than %f.\n", lb_rho);
  3395.         (void) fflush(logFile);
  3396.         break;
  3397.  
  3398. //______________________________________________________________________________
  3399.  
  3400.     case LS_search_freq:
  3401.         (void) sscanf(line, "%*s " FDFMT, &search_freq);
  3402.         pr(logFile, "Local search will be performed with frequency %f.\n", search_freq);
  3403.         (void) fflush(logFile);
  3404.         break;
  3405.  
  3406. //______________________________________________________________________________
  3407.  
  3408.     case DPF_ANALYSIS:
  3409.         /*
  3410.         ** analysis
  3411.         */
  3412.         /* _____________________________________________________________________
  3413.         **
  3414.         ** Perform Cluster analysis on results of docking,
  3415.         ** _____________________________________________________________________
  3416.         */
  3417.         if (!command_mode) {
  3418.             analysis( Nnb, atomstuff, charge, abs_charge, qsp_abs_charge, B_calcIntElec, clus_rms_tol,
  3419.                       crdpdb, ad_energy_tables, map, econf, nruns,
  3420.                       natom, nonbondlist, nconf, ntor, sHist, FN_ligand,
  3421.                       lig_center, B_symmetry_flag, tlist, type, vt, FN_rms_ref_crds,
  3422.                       torsFreeEnergy, B_write_all_clusmem, ligand_is_inhibitor,
  3423.                       outlev,
  3424.                       ignore_inter, B_include_1_4_interactions, scale_1_4,
  3425.                       parameterArray, unbound_internal_FE,
  3426.                       info, B_use_non_bond_cutoff, B_have_flexible_residues,
  3427.                       B_rms_atoms_ligand_only);
  3428.  
  3429.             (void) fflush(logFile);
  3430.         } else {
  3431.             (void)fprintf(logFile, "NOTE: Command mode has been set, so cluster analysis cannot be performed.\n\n");
  3432.         }
  3433.         break;
  3434.  
  3435. //______________________________________________________________________________
  3436.  
  3437.     case DPF_TORSDOF:
  3438.         /*
  3439.         ** torsdof %d %f
  3440.         */
  3441.         retval = sscanf( line, "%*s %d " FDFMT, &ntorsdof, &torsdoffac );
  3442.         if (retval == 2) {
  3443.             pr( logFile, "WARNING:  The torsional DOF coefficient is now read in from the parameter file; the value specified here (%.4lf) will be ignored.\n\n", (double)torsdoffac);
  3444.         }
  3445.         pr( logFile, "Number of torsional degrees of freedom = %d\n", ntorsdof);
  3446.         pr( logFile, "Free energy coefficient for torsional degrees of freedom = %.4f", AD4.coeff_tors);
  3447.         if (parameter_library_found == 1) {
  3448.             pr( logFile, " as specified in parameter library \"%s\".\n\n", FN_parameter_library );
  3449.         } else {
  3450.             pr( logFile, ", the factory default value.\n\n");
  3451.         }
  3452.  
  3453.         torsFreeEnergy = (Real)ntorsdof * AD4.coeff_tors;
  3454.  
  3455.         pr( logFile, "Estimated loss of torsional free energy upon binding = %+.4f kcal/mol\n\n", torsFreeEnergy);
  3456.         (void) fflush(logFile);
  3457.         break;
  3458.  
  3459. //______________________________________________________________________________
  3460.  
  3461.     case DPF_INVESTIGATE:
  3462.         /*
  3463.         ** Bin energies by RMSD from reference structure
  3464.         **
  3465.         ** investigate 100000 1000000 100
  3466.         */
  3467.         (void) sscanf( line, "%*s %d %d %d", &OutputEveryNTests, &maxTests, &NumLocalTests );
  3468.         (void) fprintf( logFile, "OutputEveryNTests= %d\n", OutputEveryNTests);
  3469.         (void) fprintf( logFile, "maxTests= %d\n", maxTests );
  3470.         (void) fprintf( logFile, "NumLocalTests= %d\n\n", NumLocalTests );
  3471.         (void) investigate( Nnb, charge, abs_charge, qsp_abs_charge, B_calcIntElec,
  3472.                             crd, crdpdb, ad_energy_tables,
  3473.                             maxTests,
  3474.                             map, natom, nonbondlist, ntor,
  3475.                             outlev, tlist, type, vt, B_isGaussTorCon, US_torProfile,
  3476.                             B_isTorConstrained, B_ShowTorE, US_TorE,
  3477.                             F_TorConRange, N_con, B_symmetry_flag, FN_rms_ref_crds,
  3478.                             OutputEveryNTests, NumLocalTests, trnStep0, torStep0,
  3479.                             ignore_inter,
  3480.                             B_include_1_4_interactions, scale_1_4,
  3481.                             parameterArray, unbound_internal_FE,
  3482.                             info, B_use_non_bond_cutoff, B_have_flexible_residues );
  3483.  
  3484.         (void) fflush(logFile);
  3485.         break;
  3486.  
  3487. //______________________________________________________________________________
  3488.  
  3489.     case DPF_LIG_NOT_INHIB:
  3490.         /*
  3491.         ** ligand_is_not_inhibitor
  3492.         */
  3493.         ligand_is_inhibitor = 0;
  3494.         pr( logFile, "\nThis ligand is not an inhibitor, so dissociation constants (Kd) will be calculated, not inhibition constants (Ki).\n\n" );
  3495.         (void) fflush(logFile);
  3496.         break;
  3497.  
  3498. /*____________________________________________________________________________*/
  3499.  
  3500.     case DPF_UNBOUND:
  3501.         /*
  3502.          * unbound 0.0
  3503.          */
  3504.         (void) sscanf( line, "%*s " FDFMT, &unbound_internal_FE );
  3505.         pr(logFile, "The internal energy of the unbound state was set to %+.3lf kcal/mol\n\n", unbound_internal_FE);
  3506.         (void) fflush(logFile);
  3507.         break;
  3508.  
  3509. /*____________________________________________________________________________*/
  3510.  
  3511.     case DPF_COMPUTE_UNBOUND_EXTENDED:
  3512.         /*
  3513.          *  compute_unbound_extended
  3514.          */
  3515.         if (ntor > 0) {
  3516.  
  3517.             pr(logFile, "Computing the energy of the unbound state of the ligand,\ngiven the torsion tree defined in the ligand file.\n\n");
  3518.             (void) fflush( logFile );
  3519.  
  3520.             // The initial goal is to obtain an extended conformation of the ligand.
  3521.  
  3522.             // Step 0 // {
  3523.             //
  3524.             // Set termination criteria for unbound calculations
  3525.             // -------------------------------------------------
  3526.             //
  3527.             // Set the maximum number of energy evaluations for finding the unbound conformation
  3528.             // if num_evals is less than this, then a shorter unbound docking will be performed
  3529.             max_evals_unbound = 1000000; // 1 million
  3530.             num_evals_unbound = num_evals > max_evals_unbound ?  max_evals_unbound :  num_evals;
  3531.             // end of Step 0 // }
  3532.  
  3533.             // Step 1 // {
  3534.             //
  3535.             // Run a hybrid global-local search using the unbound energy tables (set to be repulsive-only)
  3536.             // -------------------------------------------------------------------------------------------
  3537.             //
  3538.             //  *  Turn off the use of the non-bond cutoff
  3539.             //  *  Turn off internal electrostatics
  3540.             //  *  Turn off intermolecular energy calculations
  3541.             // TODO Do not translate or rotate the ligand in unbound searches
  3542.             //
  3543.             /*
  3544.              *  Genetic Algorithm-Local search,  a.k.a.
  3545.              *  Lamarckian Genetic Algorithm
  3546.              */
  3547.             if ((GlobalSearchMethod==NULL)||(LocalSearchMethod==NULL)) {
  3548.                 prStr(error_message, "%s:  ERROR:  You must use \"set_ga\" to allocate both Global Optimization object AND Local Optimization object.\n", programname);
  3549.                 stop(error_message);
  3550.                 exit(-1);
  3551.             } else if (!B_found_elecmap) {
  3552.                 prStr(error_message, "%s:  ERROR:  no \"elecmap\" command has been specified!\n", programname);
  3553.                 stop(error_message);
  3554.                 exit(-1);
  3555.             } else if (!B_found_desolvmap) {
  3556.                 prStr(error_message, "%s:  ERROR:  no \"desolvmap\" command has been specified!\n", programname);
  3557.                 stop(error_message);
  3558.                 exit(-1);
  3559.             }
  3560.             //
  3561.             // Do not use a non-bond cutoff, this helps to produce the "most" extended conformation
  3562.             // especially with long inhibitors
  3563.             B_use_non_bond_cutoff = FALSE;
  3564.             //
  3565.             // Save the current value of B_calcIntElec, so we can restore it later.
  3566.             B_calcIntElec_saved = B_calcIntElec;
  3567.             //
  3568.             // Set the calculation of internal electrostatics to FALSE:
  3569.             // B_calcIntElec = FALSE;
  3570.             //
  3571.             // Assume the unbound state of the receptor is the same as the input coordinates from the
  3572.             // flexible residues file.  This means we must not change the rotatable bonds in the
  3573.             // flexible residues of the receptor during the unbound extended search.
  3574.             // We can turn off rotation of the flexres by setting ntor to ntor_ligand.
  3575.             // Save the current value of "ntor" in the "sInit" state variable, set it to number of torsions
  3576.             // in the ligand for the unbound extended search, then restore it afterwards.
  3577.             saved_sInit_ntor = sInit.ntor;
  3578.             sInit.ntor = ntor_ligand;
  3579.             //
  3580.             // Use the repulsive unbound energy tables, "unbound_energy_tables",
  3581.             // to drive the molecule into an extended conformation
  3582.             evaluate.setup( crd, charge, abs_charge, qsp_abs_charge, type, natom, map,
  3583.                             elec, emap, nonbondlist, unbound_energy_tables, Nnb,
  3584.                             B_calcIntElec, B_isGaussTorCon, B_isTorConstrained,
  3585.                             B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, crdreo, sInit, ligand,
  3586.                             ignore_inter,
  3587.                             B_include_1_4_interactions, scale_1_4,
  3588.                             parameterArray, unbound_internal_FE, info, B_use_non_bond_cutoff, B_have_flexible_residues);
  3589.             //
  3590.             // Turn off computing the intermolecular energy, we will only consider the intramolecular energy
  3591.             // to determine the unbound state of the flexible molecule:
  3592.             evaluate.compute_intermol_energy(FALSE);
  3593.             //
  3594.             (void) fprintf( logFile, "\n\n\tBEGINNING COMPUTATION OF UNBOUND EXTENDED STATE USING LGA\n");
  3595.             (void) fprintf( logFile,     "\t_________________________________________________________\n\n\n");
  3596.             (void) fflush( logFile );
  3597.             //
  3598.             pr(logFile, "Date:\t");
  3599.             printdate( logFile, 2 );
  3600.             (void) fflush( logFile );
  3601.             gaStart = times( &tms_gaStart );
  3602.             //  Can get rid of the following line
  3603.             ((Genetic_Algorithm *)GlobalSearchMethod)->initialize(pop_size, 7+sInit.ntor);
  3604.             //
  3605.             // Start Lamarckian GA run searching only torsions -- Unbound simulation
  3606.             // sUnbound_ext = call_glss_tors( GlobalSearchMethod, LocalSearchMethod,
  3607.             sUnbound_ext = call_glss( GlobalSearchMethod, LocalSearchMethod,
  3608.                                       sInit,
  3609.                                       num_evals_unbound, pop_size,
  3610.                                       outlev,
  3611.                                       outputEveryNgens, &ligand,
  3612.                                       // B_RandomDihe0, // use this line with call_glss_tors()
  3613.                                       B_RandomTran0, B_RandomQuat0, B_RandomDihe0,
  3614.                                       info, FN_pop_file );
  3615.             // State of best individual at end of GA-LS run, sUnbound_ext, is returned.
  3616.             // Finished Lamarckian GA run
  3617.             gaEnd = times( &tms_gaEnd );
  3618.             pr( logFile, "\nFinished Lamarckian Genetic Algorithm (LGA), time taken:\n");
  3619.             timesyshms( gaEnd - gaStart, &tms_gaStart, &tms_gaEnd );
  3620.             pr( logFile, "\n");
  3621.             printdate( logFile, 1 );
  3622.             (void) fflush( logFile );
  3623.             pr(logFile, "\nTotal number of Energy Evaluations: %lu\n", evaluate.evals() );
  3624.             pr(logFile, "Total number of Generations:        %u\n", ((Genetic_Algorithm *)GlobalSearchMethod)->num_generations());
  3625.             // end of Step 1 // }
  3626.  
  3627.             // Step 2 // {
  3628.             //
  3629.             // Do a short local search using the standard internal energy tables
  3630.             // -----------------------------------------------------------------
  3631.             //
  3632.             // turn on internal electrostatics
  3633.             // but keep intermolecular energy calculations off
  3634.             //
  3635.             // Turn on calculation of internal electrostatics:
  3636.             //// B_calcIntElec = TRUE;
  3637.             //
  3638.             // Use the standard AutoDock energy tables to compute the internal energy
  3639.             // Use this value to set unbound_internal_FE
  3640.             //// evaluate.setup( crd, charge, abs_charge, qsp_abs_charge, type, natom, map,
  3641.                             //// elec, emap, nonbondlist, ad_energy_tables, Nnb,
  3642.                             //// B_calcIntElec, B_isGaussTorCon, B_isTorConstrained,
  3643.                             //// B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, crdreo, sInit, ligand,
  3644.                             //// ignore_inter,
  3645.                             //// B_include_1_4_interactions, scale_1_4,
  3646.                             //// parameterArray, unbound_internal_FE, info, B_use_non_bond_cutoff, B_have_flexible_residues);
  3647.             //
  3648.             // --- Start Local Search ---
  3649.             //// pr( logFile, "\nPerforming local search using standard AutoDock scoring function\n" );
  3650.             //// pr( logFile, "\nUsing UnboundLocalSearchMethod = new Solis_Wets1(7+sInit.ntor, 300, 4, 4, 1., 0.01, 2., 0.5, 1.);\n\n" );
  3651.             // Create a local search object
  3652.             // * Use an initial rho value of 0.1 (default is set in DPF by "sw_rho 1.0")
  3653.             //   to ensure smaller, 'more local' steps.
  3654.             // * Use a search frequency of 1.0 (default is set in DPF by "ls_search_freq 0.06")
  3655.             //// unsigned int ls_pop_size = 150;
  3656.             // max_its = 300
  3657.             // max_succ = 4
  3658.             // max_fail = 4
  3659.             // rho = 1.
  3660.             // lb_rho = 0.01
  3661.             // expansion = 2.
  3662.             // contraction = 0.5
  3663.             // search_freq = 1.
  3664.             //// UnboundLocalSearchMethod = new Solis_Wets1(7+sInit.ntor, 300, 4, 4, 1., 0.01, 2., 0.5, 1.);
  3665.             // Perform a local search, using the standard AutoDock 4 scoring function
  3666.             //// sUnbound_ls = call_ls( UnboundLocalSearchMethod, sUnbound_ext, ls_pop_size, &ligand );
  3667.             //// // sUnbound_ext = sUnbound_ls; // if you want to update sUnbound_ext to be sUnbound_ls...
  3668.             // --- Finished Local Search ---
  3669.             // end of Step 2 // }
  3670.  
  3671.             // Step 3 // {
  3672.             //
  3673.             // Restore the AutoDock 4 force field for docking
  3674.             // ----------------------------------------------
  3675.             //
  3676.             // Remember to turn on the use of the non-bond cutoff
  3677.             B_use_non_bond_cutoff = TRUE;
  3678.             //
  3679.             // Restore the setting for calculation of internal electrostatics to the saved value:
  3680.             B_calcIntElec = B_calcIntElec_saved;
  3681.             //
  3682.             // Restore the number of torsions in the state variables "sInit" and "sUnbound_ext"
  3683.             sInit.ntor = saved_sInit_ntor;
  3684.             sUnbound_ext.ntor = saved_sInit_ntor;
  3685.             //
  3686.             // Use the standard AutoDock energy tables to compute the internal energy
  3687.             // Use this value to set unbound_internal_FE
  3688.             evaluate.setup( crd, charge, abs_charge, qsp_abs_charge, type, natom, map,
  3689.                             elec, emap, nonbondlist, ad_energy_tables, Nnb,
  3690.                             B_calcIntElec, B_isGaussTorCon, B_isTorConstrained,
  3691.                             B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, crdreo, sInit, ligand,
  3692.                             ignore_inter,
  3693.                             B_include_1_4_interactions, scale_1_4,
  3694.                             parameterArray, unbound_internal_FE, info, B_use_non_bond_cutoff, B_have_flexible_residues);
  3695.             // end of Step 3 // }
  3696.  
  3697.             // Step 4 // {
  3698.             //
  3699.             // Compute the energy of the unbound extended state
  3700.             // ------------------------------------------------
  3701.             //
  3702.             // Convert from unbound state to unbound coordinates
  3703.             cnv_state_to_coords( sUnbound_ext, vt, tlist, sUnbound_ext.ntor, crdpdb, crd, natom );
  3704.             //
  3705.             // Calculate the unbound internal energy using the standard AutoDock energy function
  3706.             (void) eintcalPrint(nonbondlist, ad_energy_tables, crd, Nnb, B_calcIntElec, B_include_1_4_interactions, scale_1_4, qsp_abs_charge, parameterArray, B_use_non_bond_cutoff, B_have_flexible_residues);
  3707.             //
  3708.             // eintcal() and eintcalPrint() set the values of the hideously-global nb_group_energy[]
  3709.             unbound_ext_internal_FE = nb_group_energy[INTRA_LIGAND] + nb_group_energy[INTRA_RECEPTOR];
  3710.             //
  3711.             pr(logFile, "\n\nThe internal energy of the unbound extended state was computed to be %+.3lf kcal/mol\n\n", unbound_ext_internal_FE);
  3712.             // end of Step 4 // }
  3713.  
  3714.             // Step 5 // {
  3715.             //
  3716.             // Decide whether to use extended or AutoDock state for unbound state
  3717.             // ------------------------------------------------------------------
  3718.             //
  3719.             if (unbound_ext_internal_FE > 0.0) {
  3720.                 // Unbound extended state has an internal energy that is positive
  3721.  
  3722.                 // Step 5.1 // {
  3723.                 //
  3724.                 // Repeat Step 1 with the standard AutoDock internal energy potentials
  3725.                 //
  3726.                 // Run a hybrid global-local search using the autodock energy tables
  3727.                 // -----------------------------------------------------------------
  3728.                 //
  3729.                 //  *  Turn off the use of the non-bond cutoff
  3730.                 //  *  Turn off internal electrostatics
  3731.                 //  *  Turn off intermolecular energy calculations
  3732.                 // TODO Do not translate or rotate the ligand in unbound searches
  3733.                 //
  3734.                 /*
  3735.                  *  Genetic Algorithm-Local search,  a.k.a.
  3736.                  *  Lamarckian Genetic Algorithm
  3737.                  */
  3738.                 (void) fprintf( logFile, "\n\n\tBEGINNING COMPUTATION OF UNBOUND AUTODOCK STATE USING LGA\n");
  3739.                 (void) fprintf( logFile,     "\t_________________________________________________________\n\n\n");
  3740.                 (void) fflush( logFile );
  3741.                 //
  3742.                 pr(logFile, "Date:\t");
  3743.                 printdate( logFile, 2 );
  3744.                 (void) fflush( logFile );
  3745.                 gaStart = times( &tms_gaStart );
  3746.                 //  Can get rid of the following line
  3747.                 ((Genetic_Algorithm *)GlobalSearchMethod)->initialize(pop_size, 7+sInit.ntor);
  3748.                 //
  3749.                 // Start Lamarckian GA run searching only torsions -- Unbound simulation
  3750.                 // sUnbound_ad = call_glss_tors( GlobalSearchMethod, LocalSearchMethod,
  3751.                 sUnbound_ad = call_glss( GlobalSearchMethod, LocalSearchMethod,
  3752.                                          sInit,
  3753.                                          num_evals_unbound, pop_size,
  3754.                                          outlev,
  3755.                                          outputEveryNgens, &ligand,
  3756.                                          B_RandomTran0, B_RandomQuat0, B_RandomDihe0,
  3757.                                          info, FN_pop_file );
  3758.                 // State of best individual at end of GA-LS run, sUnbound_ad, is returned.
  3759.                 // Finished Lamarckian GA run
  3760.                 gaEnd = times( &tms_gaEnd );
  3761.                 pr( logFile, "\nFinished Lamarckian Genetic Algorithm (LGA), time taken:\n");
  3762.                 timesyshms( gaEnd - gaStart, &tms_gaStart, &tms_gaEnd );
  3763.                 pr( logFile, "\n");
  3764.                 printdate( logFile, 1 );
  3765.                 (void) fflush( logFile );
  3766.                 pr(logFile, "\nTotal number of Energy Evaluations: %lu\n", evaluate.evals() );
  3767.                 pr(logFile, "Total number of Generations:        %u\n", ((Genetic_Algorithm *)GlobalSearchMethod)->num_generations());
  3768.                 // Restore the number of torsions in the state variable "sUnbound_ad"
  3769.                 sUnbound_ad.ntor = saved_sInit_ntor;
  3770.                 // end of Step 5.1 // }
  3771.  
  3772.                 // Step 5.2 // {
  3773.                 //
  3774.                 // Compute the energy of the unbound AutoDock state
  3775.                 // ------------------------------------------------
  3776.                 //
  3777.                 // Convert from unbound state to unbound coordinates
  3778.                 cnv_state_to_coords( sUnbound_ad, vt, tlist, sUnbound_ad.ntor, crdpdb, crd, natom );
  3779.                 //
  3780.                 // Calculate the unbound internal energy using the standard AutoDock energy function
  3781.                 (void) eintcalPrint(nonbondlist, ad_energy_tables, crd, Nnb, B_calcIntElec, B_include_1_4_interactions, scale_1_4, qsp_abs_charge, parameterArray, B_use_non_bond_cutoff, B_have_flexible_residues);
  3782.                 //
  3783.                 // eintcal() and eintcalPrint() set the values of the hideously-global nb_group_energy[]
  3784.                 unbound_ad_internal_FE = nb_group_energy[INTRA_LIGAND] + nb_group_energy[INTRA_RECEPTOR];
  3785.                 //
  3786.                 pr(logFile, "\n\nThe internal energy of the unbound AutoDock state was computed to be %+.3lf kcal/mol\n\n", unbound_ad_internal_FE);
  3787.                 // end of Step 5.2 // }
  3788.  
  3789.                 if (unbound_ad_internal_FE < unbound_ext_internal_FE) {
  3790.                     pr(logFile, "NOTE:   The AutoDock internal energy of the \"extended\" state was higher\nNOTE:   than that of the state obtained by searching using the AutoDock internal\nNOTE:   energy potentials.\nNOTE:   The unbound state was set to the AutoDock optimum state, not the \"extended\" state.\n\n");
  3791.                     unbound_internal_FE = unbound_ad_internal_FE;
  3792.                     sUnbound = sUnbound_ad;
  3793.                 } else {
  3794.                     pr(logFile, "NOTE:   Although the AutoDock internal energy of the \"extended\" state was positive, it was lower\nNOTE:   than that of the state obtained by searching using the AutoDock internal\nNOTE:   energy potentials.\nNOTE:   The unbound state was set to the \"extended\" state.\n\n");
  3795.                     unbound_internal_FE = unbound_ext_internal_FE;
  3796.                     sUnbound = sUnbound_ext;
  3797.                 }
  3798.             } else {
  3799.                 // Unbound extended state has an internal energy that is negative
  3800.                 unbound_internal_FE = unbound_ext_internal_FE;
  3801.                 sUnbound = sUnbound_ext;
  3802.                 pr(logFile, "NOTE:   The AutoDock internal energy of the \"extended\" state was negative.\n\nNOTE:   The unbound state was set to the \"extended\" state.\n\n");
  3803.             }
  3804.             //
  3805.             pr(logFile, "\n\nThe internal energy of the unbound state was set to %+.3lf kcal/mol\n\n", unbound_internal_FE);
  3806.             // end of Step 5 // }
  3807.  
  3808.             // Step 6 // {
  3809.             //
  3810.             // Convert from unbound state to unbound coordinates
  3811.             cnv_state_to_coords( sUnbound, vt, tlist, sUnbound.ntor, crdpdb, crd, natom );
  3812.             // end of Step 6 // }
  3813.  
  3814.             // Step 7 // {
  3815.             //
  3816.             // Output the coordinates of the unbound state
  3817.             pr( logFile, "\n\n\tFINAL UNBOUND STATE\n" );
  3818.             pr( logFile,     "\t___________________\n\n\n" );
  3819.             //
  3820.             writePDBQT( -1, seed,  FN_ligand, dock_param_fn, lig_center,
  3821.                         sUnbound, ntor, &eintra, &einter, natom, atomstuff,
  3822.                         crd, emap, elec,
  3823.                         charge, abs_charge, qsp_abs_charge,
  3824.                         ligand_is_inhibitor,
  3825.                         torsFreeEnergy,
  3826.                         vt, tlist, crdpdb, nonbondlist,
  3827.                         ad_energy_tables,
  3828.                         type, Nnb, B_calcIntElec,
  3829.                         map,
  3830.                         outlev,
  3831.                         ignore_inter,
  3832.                         B_include_1_4_interactions, scale_1_4, parameterArray, unbound_internal_FE,
  3833.                         info, UNBOUND, PDBQT_record, B_use_non_bond_cutoff, B_have_flexible_residues);
  3834.             // end of Step 7 // }
  3835.  
  3836.             // Step 8 // {
  3837.             //
  3838.             // Remember to reset the energy evaluator back to computing the intermolecular energy between
  3839.             // the flexible and the rigid molecules.
  3840.             evaluate.compute_intermol_energy(TRUE);
  3841.             // end of Step 8 // }
  3842.  
  3843.         } else {
  3844.             pr(logFile, "NOTE:  AutoDock cannot compute the energy of the unbound state, since the ligand is rigid.\n\n");
  3845.             pr(logFile, "NOTE:  Use the \"unbound\" command to set the energy of the unbound state, if known from a previous calculation where the ligand was treated as flexible.\n\n");
  3846.             unbound_internal_FE = 0.0L;
  3847.             pr(logFile, "\n\nThe internal energy of the unbound state was set to %+.3lf kcal/mol\n\n", unbound_internal_FE);
  3848.         }
  3849.  
  3850.         pr( logFile, UnderLine );
  3851.         (void) fflush(logFile);
  3852.         break;
  3853.  
  3854. /*____________________________________________________________________________*/
  3855.  
  3856.     case DPF_EPDB:
  3857.         /*
  3858.          *  epdb
  3859.          *
  3860.          *  Computes the energy of the ligand specified by the "move lig.pdbqt" command.
  3861.          *  Return the energy of the Small Molecule.
  3862.          *  FN_ligand must be in   PDBQT-format;
  3863.          *  flag can be:-
  3864.          *  0 = NEW, or   PDBQT-71, and
  3865.          *  1 = OLD, or   PDBQT-55 (old PDBq format).
  3866.          */
  3867.         outside = FALSE;
  3868.         atoms_outside = FALSE;
  3869.         eintra = 0.0L;
  3870.         einter = 0.0L;
  3871.         etotal = 0.0L;
  3872.  
  3873.         pr(logFile, "WARNING This command, \"epdb\", currently computes the energy of the ligand specified by the \"move lig.pdbqt\" command.\n");
  3874.         retval = sscanf(line, "%*s %s", dummy_FN_ligand);
  3875.         if (retval >= 1) {
  3876.             pr(logFile, "WARNING  -- it will not read in the PDBQT file specified on the \"epdb\" command line.\n");
  3877.         }
  3878.  
  3879.         /*
  3880.         (void) sscanf(line, "%*s %s", FN_ligand);
  3881.         pr(logFile, "epdb %s\n\n", FN_ligand);
  3882.         natom = 0;
  3883.         print_1_4_message(logFile, B_include_1_4_interactions, scale_1_4);
  3884.         //
  3885.         ligand = readPDBQT( line,
  3886.                             num_atom_types,
  3887.                             &natom,
  3888.                             crdpdb, crdreo, charge, &B_haveCharges,
  3889.                             type, bond_index,
  3890.                             pdbaname, FN_ligand, FN_flexres, B_have_flexible_residues, atomstuff, Htype,
  3891.                             &B_constrain_dist, &atomC1, &atomC2,
  3892.                             &sqlower, &squpper,
  3893.                             &ntor1, &ntor, &ntor_ligand,
  3894.                             tlist, vt,
  3895.                             &Nnb, nonbondlist,
  3896.                             jobStart, tms_jobStart, hostnm, &ntorsdof, outlev,
  3897.                             ignore_inter,
  3898.                             B_include_1_4_interactions,
  3899.                             atoms, PDBQT_record );
  3900.  
  3901.         //
  3902.         // pre-calculate some values we will need later in computing the desolvation energy
  3903.         //
  3904.         for (i=0; i<natom; i++) {
  3905.             abs_charge[i] = fabs(charge[i]);
  3906.             qsp_abs_charge[i] = qsolpar * abs_charge[i];
  3907.         }
  3908.         */
  3909.  
  3910.         // to restore the original coordinates, we must temporarily undo the centering transformation
  3911.         for ( i=0; i<true_ligand_atoms; i++ ) {
  3912.             for (xyz = 0;  xyz < SPACE;  xyz++) {
  3913.                 crdpdb[i][xyz] += lig_center[xyz];
  3914.             }
  3915.         }
  3916.         // determine if any atoms are outside the grid box
  3917.         atoms_outside = FALSE;
  3918.         for (i=0; i<natom; i++) {
  3919.             outside = is_out_grid_info(crdpdb[i][X], crdpdb[i][Y], crdpdb[i][Z]);
  3920.             if (outside) {
  3921.                 atoms_outside = TRUE;
  3922.                 (void) sprintf( message, "%s: WARNING: Atom %d (%.3f, %.3f, %.3f) is outside the grid!\n", programname, i+1, crdpdb[i][X], crdpdb[i][Y], crdpdb[i][Z] );
  3923.                 print_2x( logFile, stderr, message );
  3924.                 outside = FALSE; /* Reset outside */
  3925.             }
  3926.         }
  3927.         pr(logFile, "Number of \"true\" ligand atoms:  %d\n", true_ligand_atoms);
  3928.         //
  3929.         for (i=0;i<natom;i++) {
  3930.             if (ignore_inter[i] == 1) {
  3931.                 pr(logFile, "Special Boundary Conditions:\n");
  3932.                 pr(logFile, "____________________________\n\n");
  3933.                 pr(logFile, "AutoDock will ignore the following atoms in the input PDBQT file \nin intermolecular energy calculations:\n");
  3934.                 pr(logFile, "\n(This is because these residue atoms are at the boundary between \nflexible and rigid, and since they cannot move \nthey will not affect the total energy.)\n\n");
  3935.                 break;
  3936.             }
  3937.         }
  3938.         for (i=0;i<natom;i++) {
  3939.             if (ignore_inter[i] == 1) {
  3940.                 pr(logFile, "Atom number %d:  %s\n", i+1, atomstuff[i] );
  3941.             }
  3942.         }
  3943.         pr(logFile, "\n");
  3944.  
  3945.         sInit.ntor = ligand.S.ntor;
  3946.  
  3947.         // save any currently-computed unbound internal FE
  3948.         unbound_internal_FE_saved = unbound_internal_FE;
  3949.  
  3950.         // Initialise to zero, since we may call "epdb" more than once in a single DPF
  3951.         // Set the unbound free energy -- assume it is zero, since this is a "single-point energy calculation"
  3952.         unbound_internal_FE = 0.0;
  3953.  
  3954.         // Calculate the internal energy
  3955.         if (ntor > 0) {
  3956.             (void) eintcalPrint(nonbondlist, ad_energy_tables, crdpdb, Nnb, B_calcIntElec, B_include_1_4_interactions, scale_1_4, qsp_abs_charge, parameterArray, B_use_non_bond_cutoff, B_have_flexible_residues);
  3957.         }
  3958.  
  3959.         // calculate the energy breakdown for the input coordinates, "crdpdb"
  3960.         eb = calculateEnergies( natom, ntor, unbound_internal_FE, torsFreeEnergy, B_have_flexible_residues,
  3961.                                 crdpdb, charge, abs_charge, type, map, info, outside,
  3962.                                 ignore_inter, elec, emap, &elec_total, &emap_total,
  3963.                                 nonbondlist, ad_energy_tables, Nnb, B_calcIntElec,
  3964.                                 B_include_1_4_interactions, scale_1_4, qsp_abs_charge, parameterArray, B_use_non_bond_cutoff);
  3965.  
  3966.         pr(logFile, "\n\n\t\tIntermolecular Energy Analysis\n");
  3967.         pr(logFile,     "\t\t==============================\n\n");
  3968.         pr(logFile, "Atom  vdW+Hb+Elec  vdW+Hbond  Electrosta  Partial          Coordinates         \n");
  3969.         pr(logFile, "Type    Energy      Energy    tic Energy  Charge      x         y         z    \n");
  3970.         pr(logFile, "____  __________  __________  __________  _______  ________  ________  ________\n");
  3971.         //          "1234  0123456789  0123456789  0123456789  1234567  12345678  12345678  12345678"
  3972.  
  3973.         charge_total = 0.;
  3974.         etot = 0.;
  3975.         for (i = 0;  i < natom;  i++) {
  3976.             etot = emap[i] + elec[i];
  3977.             pr(logFile, "%4d  %10.2f  %10.2f  %10.2f  %7.3f  %8.4f  %8.4f  %8.4f\n", (type[i]+1), etot, emap[i], elec[i], charge[i], crdpdb[i][X], crdpdb[i][Y], crdpdb[i][Z]);
  3978.             charge_total += charge[i];
  3979.         } /*i*/
  3980.         pr(logFile, "      __________  __________  __________  _______\n");
  3981.         pr(logFile, "Total %10.2lf  %10.2lf  %10.2lf  %7.3lf\n",        (double)(emap_total + elec_total), (double)emap_total, (double)elec_total, (double)charge_total);
  3982.         pr(logFile, "      __________  __________  __________  _______\n");
  3983.         pr(logFile, "      vdW+Hb+Elec  vdW+Hbond  Electrosta  Partial\n");
  3984.         pr(logFile, "        Energy      Energy    tic Energy  Charge\n\n");
  3985.  
  3986.         pr(logFile, "Total Intermolecular Interaction Energy   = %+.3lf kcal/mol\n", (double)(eb.e_inter + eb.e_intra) );
  3987.         pr(logFile, "Total Intermolecular vdW + Hbond Energy   = %+.3lf kcal/mol\n", (double)emap_total);
  3988.         pr(logFile, "Total Intermolecular Electrostatic Energy = %+.3lf kcal/mol\n\n\n", (double)elec_total);
  3989.  
  3990.         printEnergies( &eb, "epdb: USER    ", ligand_is_inhibitor, emap_total, elec_total, B_have_flexible_residues );
  3991.         pr(logFile, "\n");
  3992.  
  3993.         // remember to re-center the ligand
  3994.         for ( i=0; i<true_ligand_atoms; i++ ) {
  3995.             for (xyz = 0;  xyz < SPACE;  xyz++) {
  3996.                 crdpdb[i][xyz] -= lig_center[xyz];
  3997.             }
  3998.         }
  3999.  
  4000.         // restore the saved unbound internal FE
  4001.         unbound_internal_FE = unbound_internal_FE_saved;
  4002.  
  4003.         (void) fflush(logFile);
  4004.         break;
  4005.  
  4006. /*____________________________________________________________________________*/
  4007.  
  4008.     case DPF_TERMINATION:
  4009.         /*
  4010.          *  ga_termination energy 0.1
  4011.          *  ga_termination evals 25000  // the best energy did not change in this time
  4012.          *  ga_termination time 120 s
  4013.          */
  4014.         /*
  4015.         (void) sscanf( line, "%*s %d", &i );
  4016.         */
  4017.         (void) fflush(logFile);
  4018.         break;
  4019.  
  4020. /*____________________________________________________________________________*/
  4021.  
  4022.     case GA_CROSSOVER_MODE:
  4023.         /*
  4024.          * ga_crossover_mode OnePt
  4025.          * ga_crossover_mode TwoPt
  4026.          * ga_crossover_mode Uniform
  4027.          * ga_crossover_mode Arithmetic
  4028.          *
  4029.          * Xover_Mode c_mode = OnePt;  //  can be: OnePt, TwoPt, Uniform or Arithmetic
  4030.          */
  4031.         (void) sscanf( line, "%*s %s", c_mode_str );
  4032.         if (strcmp(c_mode_str, "onept") == 0) {
  4033.             c_mode = OnePt;
  4034.             pr(logFile, "One-point crossover will be used in GA and LGA searches.\n");
  4035.         } else if (strcmp(c_mode_str, "twopt") == 0) {
  4036.             c_mode = TwoPt;
  4037.             pr(logFile, "Two-point crossover will be used in GA and LGA searches.\n");
  4038.         } else if (strcmp(c_mode_str, "uniform") == 0) {
  4039.             c_mode = Uniform;
  4040.             pr(logFile, "Uniform crossover will be used in GA and LGA searches.\n");
  4041.         } else if (strcmp(c_mode_str, "arithmetic") == 0) {
  4042.             c_mode = Arithmetic;
  4043.             pr(logFile, "Arithmetic crossover will be used in GA and LGA searches.\n");
  4044.         } else {
  4045.             c_mode = Uniform; // default
  4046.         }
  4047.         (void) fflush(logFile);
  4048.         break;
  4049.  
  4050. /*____________________________________________________________________________*/
  4051.  
  4052.     case DPF_POPFILE:
  4053.         /*
  4054.          *  output_pop_file
  4055.          *
  4056.          *  Used to write out the population to a file at the end of
  4057.          *  every GA.
  4058.          */
  4059.         (void) sscanf( line, "%*s %s", FN_pop_file);
  4060.         (void) fflush(logFile);
  4061.         pr( logFile, "The population will be written to the file \"%s\" at the end of every generation.\n", FN_pop_file);
  4062.         break;
  4063.  
  4064.  /*____________________________________________________________________________*/
  4065.  
  4066.      case DPF_CONFSAMPLER:
  4067.         /*
  4068.          * confsampler
  4069.          *
  4070.          * Scan a region around conformations saved in sHist array.
  4071.          *
  4072.          */
  4073.  
  4074.         (void) sscanf( line, "%*s %s %d", confsampler_type, &confsampler_samples);
  4075.         pr( logFile, "Scanning local regions around each docked conformation.\n");
  4076.  
  4077.         if (strcmp(confsampler_type, "systematic") == 0) {
  4078.             systematic_conformation_sampler(sHist, nconf, vt, crdpdb, tlist, lig_center, natom, type, info);
  4079.         }
  4080.        
  4081.         else {
  4082.             random_conformation_sampler(sHist, nconf, confsampler_samples, vt, crdpdb, tlist, lig_center, natom, type, info);
  4083.         }
  4084.         (void) fflush(logFile);
  4085.  
  4086.         break;
  4087.  
  4088. /*____________________________________________________________________________*/
  4089.  
  4090.     case DPF_COPYRIGHT:
  4091.         /*
  4092.          * 'copyright' to show the Gnu GPL copyright
  4093.          */
  4094.         show_copyright(logFile);
  4095.         (void) fflush(logFile);
  4096.         break;
  4097.  
  4098. /*____________________________________________________________________________*/
  4099.  
  4100.     case DPF_WARRANTY:
  4101.         /*
  4102.          * 'warranty' to show the Gnu GPL warranty
  4103.          */
  4104.         show_warranty(logFile);
  4105.         (void) fflush(logFile);
  4106.         break;
  4107.  
  4108. /*_12yy_______________________________________________________________________*/
  4109.  
  4110.     case DPF_:
  4111.         /*
  4112.          *
  4113.          */
  4114.         /*
  4115.         (void) sscanf( line, "%*s %d", &i );
  4116.         (void) fflush(logFile);
  4117.         break;
  4118.         */
  4119.  
  4120. //______________________________________________________________________________
  4121.  
  4122.     default:
  4123.         /*
  4124.         **  Do nothing...
  4125.         */
  4126.         break;
  4127.  
  4128. //______________________________________________________________________________
  4129.  
  4130.     } /* switch( dpf_keyword ) */
  4131.  
  4132. } /* while PARSING-DPF parFile */
  4133.  
  4134. /* __________________________________________________________________________
  4135. **
  4136. ** Close the docking parameter file...
  4137. ** __________________________________________________________________________
  4138. */
  4139. pr( logFile, ">>> Closing the docking parameter file (DPF)...\n\n" );
  4140. //pr( logFile, UnderLine );
  4141. (void) fclose( parFile );
  4142.  
  4143.  
  4144. /* _________________________________________________________________________
  4145. **
  4146. ** If in command-mode, set the command file-pointers to standard i/o,
  4147. ** _________________________________________________________________________
  4148. */
  4149. if (command_mode) {
  4150.     status = cmdmode( natom,jobStart,tms_jobStart,
  4151.                       map, ad_energy_tables, WallEnergy, vt, tlist, ntor,
  4152.                       Nnb, nonbondlist, atomstuff, crdpdb,
  4153.                       hostnm, type, charge, abs_charge, qsp_abs_charge, B_calcIntElec,
  4154.                       atm_typ_str, torsFreeEnergy,
  4155.                       ligand_is_inhibitor,
  4156.                       ignore_inter,
  4157.                       B_include_1_4_interactions, scale_1_4,
  4158.                       parameterArray, unbound_internal_FE,
  4159.                       info, B_have_flexible_residues, B_use_non_bond_cutoff );
  4160.     exit( status );  /* "command_mode" exits here... */
  4161. }
  4162.  
  4163. //______________________________________________________________________________
  4164. /*
  4165. ** Print the time and date when the docking has finished...
  4166. */
  4167.  
  4168. pr( logFile, "This docking finished at:\t\t\t" );
  4169. printdate( logFile, 1 );
  4170. pr( logFile, "\n\n\n" );
  4171.  
  4172. success( hostnm, jobStart, tms_jobStart );
  4173.  
  4174.  if(write_stateFile){
  4175.    fprintf(stateFile,"</autodock>\n");
  4176.    (void) fclose( stateFile );
  4177.  }
  4178. (void) fclose( logFile );
  4179.  
  4180. // delete arrays
  4181. delete []nonbondlist;
  4182.  
  4183.  
  4184. //________________________________________________________________________________
  4185. /*
  4186. ** End of Boinc
  4187. */
  4188. #ifdef BOINCCOMPOUND
  4189.  boinc_fraction_done(1.);
  4190. #endif
  4191. #ifdef BOINC          
  4192.     boinc_finish(0);       /* should not return */
  4193. #endif
  4194.  
  4195. return 0;
  4196.  
  4197. } /* END OF PROGRAM */
  4198.  
  4199. #ifdef BOINC
  4200. /*  Dummy graphics API entry points.
  4201.  *  This app does not do graphics, but it still must provide these callbacks.
  4202.  */
  4203.  
  4204. void app_graphics_render(int xs, int ys, double time_of_day) {}
  4205. void app_graphics_reread_prefs(){}
  4206. void boinc_app_mouse_move(int x, int y, bool left, bool middle, bool right ){}
  4207. void boinc_app_mouse_button(int x, int y, int which, bool is_down){}
  4208. void boinc_app_key_press(int wParam, int lParam){}
  4209. void boinc_app_key_release(int wParam, int lParam){}
  4210. #endif
  4211.  
  4212. /* EOF */
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!