-
Notifications
You must be signed in to change notification settings - Fork 1
/
main.cpp
3678 lines (3264 loc) · 162 KB
/
main.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*
$Id: main.cpp,v 1.160 2021/04/04 18:10:56 mp Exp $
AutoGrid
Copyright (C) 2009 The Scripps Research Institute. All rights reserved.
AutoGrid is a Trade Mark of The Scripps Research Institute.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <sys/param.h>
#include <unistd.h> /* long sysconf(int name) */
#include <sys/types.h>
#include <ctype.h> /* tolower,isascii,isdigit */
#ifdef NOTNEEDED
#ifdef _WIN32
#include <Winsock2.h>
#include "util.h"
#endif
#endif
#include <math.h>
#include <assert.h>
#include <stdio.h>
#include <search.h>
#include <string.h>
#include <strings.h> // for bzero() on Solaris
#include <time.h>
#include <stdlib.h>
#ifndef HAVE_SYSCONF
#include "mingw_sysconf.h" // for sysconf(_SC_CLK_TCK) and possibly gethostname
#endif
#include <stddef.h>
#include <ctype.h>
#include <cassert>
/* the BOINC API header file */
#ifdef BOINC
#include "diagnostics.h"
#include "boinc_api.h"
#include "filesys.h" // boinc_fopen(), etc... */
#endif
#include "autogrid.h"
#include "autoglobal.h"
#include "autocomm.h"
#include "bondmanager.h"
#include "constants.h"
#include "distdepdiel.h"
#include "memalloc.h" // malloc_t() and calloc_t()
#include "read_parameter_library.h"
#include "threadlog.h"
#include "timesys.h"
#include "timesyshms.h"
// M Sanner 2015 add bhtree to speed up calculation using spatial hashing
#include "bhtree.h"
extern Real idct;
// round() is a C99 function and not universally available
// Required to round %.3f consistently on different platforms
#ifdef HAVE_ROUND
#define round3dp(x) ((round((x)*1000.0L))/1000.0L)
#else
#define round3dp(x) (( floor((x)*1000.0 + 0.5)) / 1000.0)
#endif
// convenience macro for plural counts in log files:
#define plural(i) ( (i)==1?"":"s" )
// convenience macro for string equality
#define streq(a,b) (strcmp((a),(b))==0)
// macros and private functions for 3-D distance work:
#define xyzxyzdist(a,b) ( hypot( ((a)[X]-(b)[X]), hypot(((a)[Y]-(b)[Y]), ((a)[Z]-(b)[Z]))))
#define vect3len(a) ( hypot( (a)[X], hypot( (a)[Y], (a)[Z])))
static double vect_sub_normalize ( double vresult[XYZ], double v1[XYZ], double v2[XYZ] );
static double vect_normalize ( double v1[XYZ] );
// distance_gt:
// shortcut return TRUE if any X,Y,Z component > limit; else returns FALSE
// and sets dc[] to components a-b and d to actual distance.
// a, b, and dc are 3-element vectors (const a,b) ; limit and d are scalars.
// NOTE dc and d are correctly set ONLY if macro returns TRUE
// CAUTION limit must be positive or zero
#define distance_gt(a, b, limit, dc, d) (\
(dc[X] = (a)[X] - (b)[X]) > (limit) || \
dc[X] < -(limit) || \
(dc[Y] = (a)[Y] - (b)[Y]) > (limit) || \
dc[Y] < -(limit) || \
(dc[Z] = (a)[Z] - (b)[Z]) > (limit) || \
dc[Z] < -(limit) , ((d = vect3len(dc))>(limit)))
// distance_le:
// shortcut return FALSE if any X,Y,Z component > limit; else returns TRUE
// and sets dc[] to components a-b and d to actual distance.
// a, b, and dc are 3-element vectors (const a,b); limit and d are scalars.
// NOTE dc and d are correctly set ONLY if macro returns TRUE
// CAUTION limit must be positive or zero
#define distance_le(a, b, limit, dc, d) (\
(dc[X] = (a)[X] - (b)[X]) <= (limit) && \
dc[X] >= -(limit) && \
(dc[Y] = (a)[Y] - (b)[Y]) <= (limit) && \
dc[Y] >= -(limit) && \
(dc[Z] = (a)[Z] - (b)[Z]) <= (limit) && \
dc[Z] >= -(limit) , ((d = vect3len(dc)) <=(limit)))
// print_error() is used with error_level where
// error_level is defined in autogrid.h
void print_error( FILE *logFile, int error_level, char *message)
// print an error or informational message to a file-pointer or
// standard error
{
char output_message[LINE_LEN];
char *tag;
switch ( error_level ) {
case FATAL_ERROR:
tag = "ERROR";
break;
case AG_ERROR:
case WARNING:
tag = "WARNING";
break;
default:
case INFORMATION:
tag = "INFORMATION";
break;
case SUGGESTION:
tag = "SUGGESTION";
break;
}
(void) snprintf( output_message, sizeof output_message,
"\n%s: %s: %s\n", programname, tag, message);
// Records all messages in the logFile.
(void) fprintf( logFile, "%s\n", output_message);
fflush(logFile);
// send only errors and fatal errors to stderr, not warnings
switch ( error_level ) {
case AG_ERROR:
case FATAL_ERROR:
(void) fprintf( stderr, "%s\n", output_message);
break;
}
// If this is a fatal error, exit now.
if (error_level == FATAL_ERROR) {
(void) fprintf( logFile, "\n%s: Unsuccessful Completion.\n", programname);
exit( EXIT_FAILURE ); // POSIX, defined in stdlib.h (usually 1)
}
}
/* fopen rewrite to either use BOINC api or normal system call */
FILE *ad_fopen(const char *path, const char *mode, FILE *logFile)
{
FILE *filep;
#ifdef BOINC
int rc;
char resolved_name[512];
rc = boinc_resolve_filename(path, resolved_name, sizeof(resolved_name));
if (rc){
fprintf(stderr, "BOINC_ERROR: cannot open filename.%s\n",path);
boinc_finish(rc); /* back to BOINC core */
}
// Then open the file with boinc_fopen() not just fopen()
filep = boinc_fopen(resolved_name, mode);
#else
filep = fopen(path,mode);
#endif
return filep;
}
static int get_rec_index(const char key[]);
// to support use_vina_potential
static int get_map_index(const char key[]);
#ifdef _OPENMP
/* M Pique */
#include <omp.h>
// MAXTHREADS is max number of hardware threads to use for computation.
#define MAXTHREADS 32
#else
#define omp_get_thread_num() (0)
#define omp_get_max_threads() (1)
#define MAXTHREADS 1
#endif
int main( int argc, char **argv )
/******************************************************************************/
/* Name: main (executable's name is "autogrid4"). */
/* Function: Calculation of interaction energy grids for Autodock. */
/* Directional H_bonds from Goodford: */
/* Distance dependent dielectric after Mehler and Solmajer. */
/* Charge-based desolvation */
/*Copyright (C) 2009 The Scripps Research Institute. All rights reserved. */
/* */
/* Authors: Garrett Matthew Morris, Ruth Huey, David S. Goodsell */
/* */
/* The Scripps Research Institute */
/* Department of Molecular Biology, MB5 */
/* 10550 North Torrey Pines Road */
/* La Jolla, CA 92037-1000. */
/* */
/* e-mail: garrett@scripps.edu */
/* rhuey@scripps.edu */
/* goodsell@scripps.edu */
/* */
/* Helpful suggestions and advice: */
/* Arthur J. Olson */
/* Bruce Duncan, Yng Chen, Michael Pique, Victoria Roberts */
/* Lindy Lindstrom */
/* */
/* Date: 07/07/04 */
/* */
/* Inputs: Control file, receptor PDBQT file, parameter file */
/* Returns: Atomic affinity, desolvation and electrostatic grid maps. */
/* Globals: NDIEL, MAX_MAPS */
/* increased from 8 to 16 6/4/2004 */
/* */
/* Modification Record */
/* Date Inits Comments */
/* 07/06/89 DSG FORTRAN implementation */
/* 07/05/92 GMM C translation */
/* 20/09/95 GMM/DSG AutoGrid3 */
/* 07/07/04 DSG/RH AutoGrid4 */
/******************************************************************************/
/* Note: 21/03/03 GMM note: ATOM_MAPS is no longer used here; was used for
* is_covalent and is_hbonder, but these are now folded into the MapObject
* and arrayed up to MAX_MAPS (currently). MAX_MAPS is always larger than
* ATOM_MAPS, so this is safe. */
{
/* use vina potential function instead of autodock4 potential function for grid calculations */
/* EXPERIMENTAL FUNCTION: NOT SUPPORTED */
const int use_vina_potential = FALSE;
/* for associative dictionary storing parameters by autogrid 'type' */
/*see atom_parameter_manager.c */
static ParameterEntry thisparm;
ParameterEntry * found_parm;
static char FN_parameter_library[MAX_CHARS]; // the AD4 parameters .dat file name
bool parameter_library_found = FALSE;
/* LIGAND:
* maximum is MAX_MAPS (really ATOM_MAPS) */
/*each type is now at most two characters plus '\0'*/
char ligand_types[MAX_MAPS][3]; /* one entry for each atom map, so two for separate_desolvation_maps */
/*array of ptrs used to parse input line*/
char * ligand_atom_types[MAX_MAPS];
/*malloc this after the number of receptor types is parsed*/
static EnergyTables et;
/* the mapObject gridmap[] holds metadata about each map.
AutoGrid map working storage is in "energy", allocated later
*/
typedef struct mapObject {
int atom_type; /*corresponds to receptor numbers????*/
int map_index;
int is_covalent;
int is_hbonder;
FILE *map_fileptr;
char map_filename[MAX_CHARS];
char type[3]; /*eg HD or OA or NA or N*/
double energy_max;
double energy_min;
double *energy; // allocated n1[X]*n1[Y]*n1[Z]*sizeof(map element [double])
double vol_probe;
double solpar_probe;
/*new 6/28*/
double Rij;
double epsij; /* already weighted by coeff_vdW */
hbond_type hbond; /*hbonding character: */
double Rij_hb;
double epsij_hb; /* already weighted by coeff_hbond */
/*per receptor type parameters, ordered as in receptor_types*/
double cA[NUM_RECEPTOR_TYPES], cB[NUM_RECEPTOR_TYPES];/*coefficients if specified in gpf*/
double nbp_r[NUM_RECEPTOR_TYPES]; /*radius of energy-well minimum*/
double nbp_eps[NUM_RECEPTOR_TYPES];/*depth of energy-well minimum*/
int xA[NUM_RECEPTOR_TYPES]; /*generally 12*/
int xB[NUM_RECEPTOR_TYPES]; /*6 for non-hbonders 10 for h-bonders*/
int hbonder[NUM_RECEPTOR_TYPES];
} MapObject;
MapObject *gridmap = NULL; /* was statically assigned MapObject gridmap[MAX_MAPS]; */
// convenience macro for indexing into energy array
// Z varies slowest (plane), then Y (rows), then X fastest (columns)
// usage: gridmap[mapnum].energy[(mapindex(ix,iy,iz)] = value;
#define mapindex2(ix,iy) ( (ix) + ((iy)*n1[X]) )
#define mapindex(ix,iy,iz) ( (ix) + n1[X] * ( (iy) + n1[Y] * (iz) ) )
// two kinds of "distance" grids, distinct from the affinity grids
static double *r_min; /* allocated full [z][y][x] for floating grid */
static int *c_count; /* allocated full [z][y][x] for constriction grid */
/*variables for RECEPTOR:*/
/*each type is now at most two characters, eg 'NA\0'*/
/*NB: these are sparse arrays, some entries are not set */
char receptor_types[NUM_RECEPTOR_TYPES][3];
/* number of different receptor atom types declared on receptor_types line in GPF */
int receptor_types_gpf_ct = 0;
int has_receptor_types_in_gpf = 0;
/* number of different receptor atom types actually found in receptor PDBQT */
int receptor_types_ct = 0;
/* array of numbers of each type */
/*NB: this is a sparse int array, some entries are 0*/
int receptor_atom_type_count[NUM_RECEPTOR_TYPES];
int lc; /* receptor pdbqt file line counter */
/*array of ptrs used to parse input line*/
char * receptor_atom_types[NUM_RECEPTOR_TYPES];
// M Sanner 2015 BHTREE
const double BH_collision_dist=1.5; // supported only with USE_BHTREE so far, but defined here anyway for simplicity
#ifdef USE_BHTREE
BHtree *bht;
BHpoint **BHat;
// MP next is strictly for debug
//#define BH_collision_dist 3.0
#define BH_cutoff_dist NBC
int *closeAtomsIndices[MAXTHREADS];
float *closeAtomsDistances[MAXTHREADS];
// M Sanner 2015 BHTREE END
#endif
/* AG_MAX_ATOMS */
/* changed these from "double" to "static" to reduce stack usage - MPique 2012 */
static double charge[AG_MAX_ATOMS];
static double vol[AG_MAX_ATOMS];
static double solpar[AG_MAX_ATOMS];
/*integers are simpler!*/
static int atom_type[AG_MAX_ATOMS];
static hbond_type hbond[AG_MAX_ATOMS];
static bool disorder[AG_MAX_ATOMS];
static char * hbtname[] = { "D0", "DS", "D1", "AS", "A1", "A2", "AD", "??" }; /* for table printouts only */
// two arrays added 2018-11 for tracking h-bond donor/acceptor neighborhoods
// TODO MP make these dynamic and allocated only if disorder_h flag is true
static int nbonds[AG_MAX_ATOMS];
static int bonded[AG_MAX_ATOMS][AG_MAX_NBONDS];
static int rexp[AG_MAX_ATOMS];
static double coord[AG_MAX_ATOMS][XYZ];
// two arrays that hold a 3-D vector for each atom (often empty):
static double rvector[AG_MAX_ATOMS][XYZ];
static double rvector2[AG_MAX_ATOMS][XYZ];
/*
(Comments by M Pique 2018-11-20 based on reading the existing code;
these document my reading of the code before some reorganization
to improve the disorder_h function):
Setting of rvector and rvector2 (both are initally empty - all zeros;
results are always attempted to be normalized):
if an atom (ib) is near an atom (ia) of class: (4 cases)
1. hbond[ia] == 2 ; D1 hydrogen bond donor (ia)
rvector[ia] is set to vector hydrogen donor TO
the other atom (ib) [i.e. ib-ia]
2. hbond[ia] == 5 ; A2 oxygen (ia)
if the oxygen has zero bonds, a warning will be issued.
rvector[ia] and rvector2[ia] will remain empty.
if the oxygen has exactly one bond it is assumed to be a carbon (atom i1)
thus the O (ia) is a Carbonyl Oxygen O=C-X
rvector[ia] is set to vector C to O [i.e., O-C: ia-i1]
and a third atom X (i2), bonded to the carbon is found.
rvector2[ia] is set to the cross product i2-i1 and O=C rvector
(C=O cross C-X gives the lone pair plane normal).
if a third atom is not found, rvector2 will be empty.
if the oxygen has exactly two bonds (i1, i2) it is Hydroxyl or Ether Oxygen X1-O-X2
The 'disorder_h' global boolean option controls the behavior.
if disorder_h and exactly one of (i1, i2) is hydrogen, atom (ib) is
set to the one that is carbon or arom_carbon (if it is).
rvector[ia] is set to the vector (ib) to O [ia-ib];
rvector2[ia] will be empty.
else [i.e., if 'not disorder_h', or neither of (i1, i2) is hydrogen,
or i1 and i2 are the same type],
rvector2[ia] is set to the vector (i1)=(i2), the lone pair plane
rvector[ia] is set to vector from the oxygen to between the lone pairs
if the oxygen has three or more bonds, only the first two are considered
and a warning will be issued.
3. hbond[ia] == 4 ; A1 directional nitrogen acceptor (ia)
if the nitrogen has zero bonds, a warning will be issued.
rvector[ia] and rvector2[ia] will remain empty.
if the nitrogen has one to three bonds
(Azide Nitrogen :N=C-X , X1-N=X2, or ...)
rvector[ia] := vector to N from mean_position(i1,i2,i3)
rvector2[ia] remains empty
if the nitrogen has four or more bonds, only the first three are considered
but no warning will be issued.
Otherwise rvector[ia] remains empty.
Use of rvector, rvector2: for each grid point, for each atom ia
if disorder_h && disorder[ia] && atom_type[ia] == hydrogen, ignore (ia) completly
For other atoms, set 'rdon','racc', 'Hramp' (for gridpoint)
based on atom (ia) class:
1. hbond[ia] == 2 ; D1 hydrogen bond donor (ia)
ia-th receptor atom = Hydrogen => receptor H-bond donor, OH or NH.
Uses dot product of rvector[ia] and rvector[closest hydrogen to grid point].
Does not use rvector2
2. hbond[ia] == 4; A1, Directional N acceptor (ia)
Uses dot product of rvector[ia] and grid point-to-ia
Does not use rvector2
3. hbond[ia] == 5 ; A2, receptor H-bond acceptor, oxygen.
if not disorder[ia] uses dot product of rvector[ia] and rvector2[ia]
if disorder[ia] "cylindrically disordered hydroxyl",
Uses dot product of rvector[ia] and grid point-to-ia
Does not use rvector2
After this point, rvector[ia] and rvector2[ia] are unused
At this grid point, for each map_index type
if gridmap[map_index].is_hbonder [.hbond>0: DS,D1,AS,A1,A2,AD]
"current map_index PROBE forms H-bonds"
sets 'rsph' based on atom ia type and map_index, range 0..1
(4 cases, with disorder a subcase of 1.)
1. if gridmap[map_index].hbond=={3,5,6} [AS or A2 or AD]
&& hbond[ia]=={1,2,6} [DS or D1 or AD]
"PROBE can be an H-BOND ACCEPTOR"
if disorder[ia]
gridmap.energy += ...[r-1.10][ hydrogen][map_index] * Hramp *
(racc + (1. - racc)*rsph)
else gridmap.energy += ...[r][atom_type[ia]][map_index] * Hramp *
(racc + (1. - racc)*rsph)
2. elseif gridmap[map_index].hbond=={4,6} [A1 or AD]
&& hbond[ia]=={1,2,6} [DS or D1 or AD]
sets gridmap{min,max,flag} from 'racc,rsph',
not depending on disorder,
does not set gridmap.energy
3. elseif gridmap[map_index].hbond=={1,2,6} [DS or D1 or AD]
&& hbond[ia]>2 [{3,4,5,6}: AS or A1 or A2 or AD]
"PROBE is H-BOND DONOR"
sets gridmap{min,max,flag} from 'rdon,rsph',
not depending on disorder,
does not set gridmap.energy
4. else
"hbonder PROBE-ia cannot form a H-bond"
sets gridmap.energy not depending on disorder.
else [if not is_hbonder]
"PROBE does not form H-bonds"
sets gridmap.energy not depending on disorder (as in case 4. above)
*/
static int hbond_12[AG_MAX_ATOMS], nhbond_12; // indices of hbond[ia]==1 or 2; count
/* receptor atom type number for autodock4 potential*/
static int hydrogen, nonHB_hydrogen, carbon, arom_carbon,
oxygen, nitrogen, nonHB_nitrogen, sulphur, nonHB_sulphur,
bromine, chlorine, fluorine, iodine;
/*canned ligand atom type number for vina_potential*/
int map_hydrogen, map_carbon, map_arom_carbon, map_oxygen, map_nitrogen;
int map_nonHB_hydrogen, map_nonHB_nitrogen, map_sulphur, map_nonHB_sulphur;
int map_chlorine, map_bromine, map_fluorine, map_iodine;
/* XYZ */
double cext[XYZ];
double cgridmax[XYZ];
double cgridmin[XYZ];
double cmax[XYZ];
double cmin[XYZ];
double csum[XYZ];
double cmean[XYZ]={0.,0.,0.}; /* receptor mean atom coords */
double center[XYZ];
double covpos[XYZ]={0.,0.,0.}; /* Cartesian-coordinate of covalent affinity well. */
int nelements[XYZ];
static int n1[XYZ]; /* nelements[i]+1: static to detect 'not set yet' */
int ne[XYZ]; /* floor (nelements[i]/2) */
/* MAX_CHARS: made static so will have strlen == 0 if not set */
static char AVS_fld_filename[MAX_CHARS];
static char floating_grid_filename[MAX_CHARS];
static char constriction_grid_filename[MAX_CHARS];
static char host_name[MAX_CHARS];
static char receptor_filename[MAX_CHARS];
static char xyz_filename[MAX_CHARS];
/* LINE_LEN */
char message[LINE_LEN];
char line[LINE_LEN];
char GPF_line[LINE_LEN];
int length = LINE_LEN;
/* NDIEL (old name MAX_DIST) for dielectric and desolvation interactions */
/* NEINT - for vdW and Hb interactions */
double energy_smooth[NEINT];
int nthreads=1; /* for OpenMP */
int iz; /* plane-counter */
char atom_name[6];
char token[LINE_LEN];
////bool warned = false;
static const char xyz[] = "xyz"; // used to print headings
static FILE *receptor_fileptr,
*AVS_fld_fileptr,
*xyz_fileptr,
*floating_grid_fileptr,
*constriction_grid_fileptr;
/*for NEW3 desolvation terms*/
double solpar_q = .01097; /*unweighted value restored 3:9:05 */
/*double solpar_q = 0.0013383; =.01097 * 0.122*/
Linear_FE_Model AD4; // set in setup_parameter_library and read_parameter_library
double q_tot = 0.0;
double diel, invdielcal=0.;//expected never used if uninitialized
double PI_halved;
double q_max = -BIG, q_min = BIG;
double r_smooth = 0.5; //NEW ON BY DEFAULT Feb2012
double spacing = 0.375; /* One quarter of a C-C bond length. */
double covhalfwidth = 1.0;
double covbarrier = 1000.0;
const double ln_half = log(0.5);
#ifndef PACKAGE_VERSION
static char * version_num = "4.2.7.x";
#else
static char * version_num = PACKAGE_VERSION;
#endif
const double factor=332.0L; /* Used to convert between calories and SI units */
/*int num_rec_types = 0;*/
float timeRemaining = 0.;
int num_atom_types = 0; /* number of ligand atom types, from "ligand_types" keyword */
int num_atom_maps = 0; /* number of ligand atom maps, from "ligand_types" keyword but larger if separate_desolvation_maps */
int num_maps = 0; /* number of "map", "elecmap", or "dsolvmap" keywords handled so far */
int elecPE = -1; /* index (num_maps value) for electrostatic map, if requested */
int dsolvPE = -1; /* index (num_maps value) for desolvation map, if requested */
bool separate_desolvation_maps = FALSE; /* write affinity and desolvation into successive maps instead of summing */
int num_distance_maps=0; // floating grid and/or constriction grid
bool floating_grid = FALSE; // Untested for a long time - M Pique 2015
bool constriction_grid = FALSE; // Added July 2020 - M Pique
double constriction_distance_cutoff = 7.5; // change with constriction_distance_cutoff GPF statement
bool dddiel = FALSE, disorder_h = FALSE;
bool map_receptor_interior = FALSE; // default with USE_BHTREE is to fast-skip over grid points within the receptor
// set the value to be reported in the atomic affinity maps for interior grid points
#define INTERIOR_VALUE 9999.
int fprintf_retval = 0;
int GPF_keyword = -1;
int indcom = 0;
int infld;
int nDone = 0;
int problem_wrt = FALSE;
//MP int xA, xB;
int outlev = LOGFORADT;
#define INIT_NUM_GRID_PTS -1
int num_grid_points_per_map = INIT_NUM_GRID_PTS;
int num_receptor_atoms=0;
static long clktck = 0;
Clock job_start;
Clock job_end;
struct tms tms_job_start;
struct tms tms_job_end;
for (int i=0; i<MAX_MAPS; i++) {
/* initialize to "" */
strcpy(ligand_types[i], "");
}
for (int i=0; i<NUM_RECEPTOR_TYPES; i++) {
/* initialize to "" */
strcpy(receptor_types[i], "");
receptor_atom_type_count[i]=0;
}
int rc;
#ifdef BOINC
int flags = 0;
flags =
BOINC_DIAG_DUMPCALLSTACKENABLED |
BOINC_DIAG_HEAPCHECKENABLED |
BOINC_DIAG_REDIRECTSTDERR |
BOINC_DIAG_REDIRECTSTDOUT ;
boinc_init_diagnostics(flags);
#ifdef BOINCCOMPOUND
BOINC_OPTIONS options;
options.main_program = false;
options.check_heartbeat = false;// monitor does check heartbeat
options.handle_trickle_ups = false;
options.handle_trickle_downs = false;
options.handle_process_control = false;
options.send_status_msgs = true;// only the worker programs (i.e. model) sends status msgs
options.direct_process_action = true;// monitor handles suspend/quit, but app/model doesn't
// Initialization of Boinc
rc = boinc_init_options(options); //return 0 for success
if( rc ){
fprintf(stderr,"BOINC_ERROR: boinc_init_options() failed \n");
exit(rc);
}
#else
// All BOINC applications must initialize the BOINC interface:
rc = boinc_init();
if (rc){
fprintf(stderr, "BOINC_ERROR: boinc_init() failed.\n");
exit(rc);
}
#endif
#endif
/*
* Get the time at the start of the run...
*/
job_start = times( &tms_job_start);
/*
* Parse the command line, open logFile and GPF file,
# report compilation option values...
*/
(void) setflags( argc, argv, version_num,
#ifdef USE_BHTREE
TRUE
#else
FALSE
#endif
,
#ifdef _OPENMP
TRUE, MAXTHREADS
#else
FALSE, 1
#endif
,
&logFile /* opened and set by setflags*/ );
/*
* Fetch clock ticks per second.
*/
if (clktck == 0) {
if ( (clktck = sysconf(_SC_CLK_TCK)) < 0) {
(void) fprintf( stderr, "\"sysconf(_SC_CLK_TCK)\" command failed in \"main.c\"\n");
(void) fprintf( logFile, "\"sysconf(_SC_CLK_TCK)\" command failed in \"main.c\"\n");
exit(EXIT_FAILURE);
} else {
idct = 1. / (float)clktck;
}
}
/* Initialize max and min coodinate bins */
for (int i = 0; i < XYZ; i++) {
cmax[i] = -BIG;
cmin[i] = BIG;
csum[i] = 0.;
}
PI_halved = PI/2.;
/*
* Initialize int receptor_atom_type_count[] array to 0
*/
for (int i=0; i<NUM_RECEPTOR_TYPES; i++) {
receptor_atom_type_count[i] = 0;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*
* Output the "AutoGrid" banner...
*/
banner( version_num, logFile);
/* report compilation options: this is mostly a duplicate of code in setflags.cpp */
(void) fprintf(logFile, " $Revision: 1.160 $\n");
(void) fprintf(logFile, "Compilation parameters: NUM_RECEPTOR_TYPES=%d NEINT=%d\n",
NUM_RECEPTOR_TYPES, NEINT);
(void) fprintf(logFile, " AG_MAX_ATOMS=%d AG_MAX_NBONDS=%d MAX_MAPS=%d NDIEL=%d MAX_ATOM_TYPES=%d\n",
AG_MAX_ATOMS, AG_MAX_NBONDS, MAX_MAPS, NDIEL, MAX_ATOM_TYPES);
fprintf(logFile," e_vdW_Hb table allows %8d entries of size %d\n",
MAX_ATOM_TYPES*MAX_ATOM_TYPES, NEINT);
/*
* Print out MAX_MAPS - maximum number of maps allowed
*/
(void) fprintf(logFile, "Maximum number of maps that can be computed = %d (defined by MAX_MAPS in \"autocomm.h\").\n", MAX_MAPS);
fprintf(logFile, " Non-bond cutoff for internal energy calculation (SOFTNBC): %.2f\n", SOFTNBC);
fprintf(logFile, " Optimize internal energy scoring (USE_8A_NBCUTOFF): ");
#ifdef USE_8A_NBCUTOFF
fprintf(logFile, " yes\n");
#else
fprintf(logFile, " no\n");
#endif
fprintf(logFile, " Faster search for nearby atoms (USE_BHTREE): ");
#ifdef USE_BHTREE
fprintf(logFile, " yes\n");
#else
fprintf(logFile, " no\n");
#endif
fprintf(logFile, " Run calculations in parallel if possible (_OPENMP): ");
#ifdef _OPENMP
fprintf(logFile, " yes\n");
fprintf(logFile, " Maximum number of parallel threads (MAXTHREADS): %d\n", MAXTHREADS);
#else
fprintf(logFile, " no\n");
#endif
/*
* Print the time and date when the file was created...
*/
(void) fprintf( logFile, "This file was created at:\t\t\t");
printdate( logFile, 1);
(void) strcpy(host_name, "unknown_host");
#ifdef HAVE_GETHOSTNAME
gethostname( host_name, sizeof host_name );
#endif
(void) fprintf( logFile, " using:\t\t\t\"%s\"\n", host_name);
(void) fprintf( logFile, "\n\n");
//______________________________________________________________________________
//
// Read in default parameters
//
setup_parameter_library(logFile, outlev, "default Unbound_Same_As_Bound", Unbound_Same_As_Bound, &AD4);
/******************************************************************************/
/* Read in the grid parameter file... */
while( fgets( GPF_line, LINE_LEN, GPF ) != NULL ) {
/******************************************************************************/
GPF_keyword = gpfparser( GPF_line);
/* This first "switch" figures out how to echo the current GPF line. */
switch( GPF_keyword ) {
case -1:
(void) fprintf( logFile, "GPF> %s", GPF_line);
print_error( logFile, WARNING, "Unrecognized keyword in grid parameter file.\n" );
continue; /* while fgets GPF_line... */
case GPF_NULL:
case GPF_COMMENT:
(void) fprintf( logFile, "GPF> %s", GPF_line);
break;
default:
(void) fprintf( logFile, "GPF> %s", GPF_line);
indcom = strindex( GPF_line, "#");
if (indcom != -1) {
GPF_line[ indcom ] = '\0'; /* Truncate str. at the comment */
}
(void) fflush( logFile);
break;
} /* first switch */
/******************************************************************************/
/* This second switch interprets the current GPF line. */
switch( GPF_keyword ) {
/******************************************************************************/
case GPF_NULL:
case GPF_COMMENT:
break;
/******************************************************************************/
case GPF_RECEPTOR:
/* read in the receptor filename */
(void) sscanf( GPF_line, "%*s %s", receptor_filename);
if ( 0==strlen(receptor_filename)) {
print_error( logFile, FATAL_ERROR, "No name specified for receptor filename");
}
(void) fprintf( logFile, "\nReceptor Input File :\t%s\n\nReceptor Atom Type Assignments:\n\n", receptor_filename);
/* try to open receptor file */
if ( (receptor_fileptr = ad_fopen(receptor_filename, "r", logFile)) == NULL ) {
(void) snprintf( message, sizeof message, "can't find or open receptor PDBQT file \"%s\".\n", receptor_filename);
print_error( logFile, FATAL_ERROR, message );
}
/* start to read in the lines of the receptor file */
num_receptor_atoms = 0;
lc = 0; /* receptor pdbqt file line counter */
while ( (fgets(line, length, receptor_fileptr)) != NULL ) {
lc++;
if (equal(line, "ATOM ") || /* Amino Acid or DNA/RNA atoms */
equal(line, "HETATM") || /* Non-standard heteroatoms */
equal(line, "CHAR")) { /* Partial Atomic Charge - not a PDB record */
/* Check that there aren't too many atoms... */
if (num_receptor_atoms >= AG_MAX_ATOMS) {
(void) sprintf( message, "Too many atoms in receptor PDBQT file %s line %d;", receptor_filename, lc );
print_error( logFile, AG_ERROR, message );
(void) sprintf( message, "-- the maximum number of atoms, AG_MAX_ATOMS, allowed is %d.", AG_MAX_ATOMS );
print_error( logFile, AG_ERROR, message );
(void) sprintf( message, "Increase the value in the \"#define AG_MAX_ATOMS %d\" line", AG_MAX_ATOMS );
print_error( logFile, SUGGESTION, message );
print_error( logFile, SUGGESTION, "in the source file \"autogrid.h\", and re-compile AutoGrid." );
(void) fflush( logFile);
// FATAL_ERROR will cause AutoGrid to exit...
print_error( logFile, FATAL_ERROR, "Sorry, AutoGrid cannot continue.");
} /* endif */
/* Check that line is long enough */
if (strlen(line) < 78) {
(void) sprintf( message,
"ATOM/HETATM line is too short in receptor PDBQT file %s line %d;",
receptor_filename, lc );
print_error( logFile, FATAL_ERROR, message );
}
(void) strncpy( atom_name, &line[12], 4);
/* atom_name is declared as an array of 6 characters,
* the PDB atom name is 4 characters (C indices 0, 1, 2 and 3)
* but let's ensure that the fifth character (C index 4)
* is a null character, which terminates the string. */
atom_name[4] = '\0';
/* Output the serial number of this atom... */
if(outlev>=LOGRECREAD)
(void) fprintf( logFile, "Atom no. %2d, \"%s\"", num_receptor_atoms + 1, atom_name);
/* Read in this receptor atom's coordinates,partial charges, and
* solvation parameters in PDBQS format... */
char field[10], field1[10], field2[10], field3[10];
(void) strncpy(field1, &line[30], 8); field[8] = '\0';
(void) strncpy(field2, &line[38], 8); field[8] = '\0';
(void) strncpy(field3, &line[46], 8); field[8] = '\0';
if ( 3 !=
sscanf(field1, "%lf", &coord[num_receptor_atoms][X]) +
sscanf(field2, "%lf", &coord[num_receptor_atoms][Y]) +
sscanf(field3, "%lf", &coord[num_receptor_atoms][Z]) ) {
(void) sprintf( message, "ATOM/HETATM line bad x,y,z in receptor PDBQT file %s line %d;", receptor_filename, lc);
print_error( logFile, FATAL_ERROR, message );
}
/* Output the coordinates of this atom... */
if(outlev>=LOGRECREAD)
(void) fprintf( logFile, " at (%.3lf, %.3lf, %.3lf), ",
coord[num_receptor_atoms][X], coord[num_receptor_atoms][Y], coord[num_receptor_atoms][Z]);
/*1:CHANGE HERE: need to set up vol and solpar*/
(void) strncpy(field, &line[70], 6); field[6] = '\0';
(void) sscanf(field, "%lf", &charge[num_receptor_atoms]);
//printf("new type is: %s\n", &line[77]);
(void) strncpy(field, &line[77], 2); field[2] = '\0';
(void) sscanf(field, "%s", thisparm.autogrid_type);
found_parm = apm_find(thisparm.autogrid_type);
if ( found_parm != NULL ) {
//(void) fprintf ( logFile, "DEBUG: found_parm->rec_index = %d, ->xs_radius= %f", found_parm->rec_index, found_parm->xs_radius);
if ( found_parm->rec_index < 0 ) {
strcpy( receptor_types[ receptor_types_ct ], found_parm->autogrid_type );
found_parm->rec_index = receptor_types_ct++;
//(void) fprintf ( logFile, "DEBUG: found_parm->rec_index => %d", found_parm->rec_index );
}
atom_type[num_receptor_atoms] = found_parm->rec_index;
solpar[num_receptor_atoms] = found_parm->solpar;
vol[num_receptor_atoms] = found_parm->vol;
hbond[num_receptor_atoms] = found_parm->hbond; /*NON=0, DS,D1, AS, A1, A2, AD */ /* N3P: added AD*/
// build list of "hbond" atoms to speed search M Pique Oct 2015
if(hbond[num_receptor_atoms]==1 || hbond[num_receptor_atoms]==2)
hbond_12[nhbond_12++] = num_receptor_atoms;
#ifdef DEBUG
printf("%d:key=%s, type=%d,solpar=%f,vol=%f\n",num_receptor_atoms,thisparm.autogrid_type, atom_type[num_receptor_atoms],solpar[num_receptor_atoms],vol[num_receptor_atoms]);
#endif
++receptor_atom_type_count[found_parm->rec_index];
} else {
char message[1000];
sprintf(message,
"\n\nreceptor file contains unknown type: '%s line %d'\nadd parameters for it to the parameter library first\n",
thisparm.autogrid_type, lc);
print_error(logFile, FATAL_ERROR, message);
}
/* if from pdbqs: convert cal/molA**3 to kcal/molA**3 */
/*solpar[num_receptor_atoms] *= 0.001;*/
q_max = max(q_max, charge[num_receptor_atoms]);
q_min = min(q_min, charge[num_receptor_atoms]);
if (atom_name[0] == ' ') {
/* truncate the first character... */
atom_name[0] = atom_name[1];
atom_name[1] = atom_name[2];
atom_name[2] = atom_name[3];
atom_name[3] = '\0';
} else if (isascii(atom_name[0])&& isdigit(atom_name[0]) &&
atom_name[1] == 'H' ) {
/* Assume this is the 'mangled' name of a hydrogen atom,
* after the atom name has been changed from 'HD21' to '1HD2'
* for example.
*
* [0-9]H\(.\)\(.\)
* 0 1 2 3
* : : : :
* V V V V
* tmp 0 1 2
* tmp
* :
* V
* 0 1 2 3
* : : : :
* V V V V
* H\(.\)\(.\)[0-9]
*/
char temp_char = atom_name[0];
atom_name[0] = atom_name[1];
atom_name[1] = atom_name[2];
atom_name[2] = atom_name[3];
atom_name[3] = temp_char;
}
/* Tell the user what you thought this atom was... */
if(outlev>=LOGRECREAD)
(void) fprintf( logFile, " was assigned atom type \"%s\" (rec_index= %d, atom_type= %d).\n", found_parm->autogrid_type, found_parm->rec_index, atom_type[num_receptor_atoms]);
/* Count the number of each atom type */
/*++receptor_atom_type_count[ atom_type[num_receptor_atoms] ];*/
/* Keep track of the extents of the receptor */
for (int i = 0; i < XYZ; i++) {
cmax[i] = max(cmax[i], coord[num_receptor_atoms][i]);
cmin[i] = min(cmin[i], coord[num_receptor_atoms][i]);
csum[i] += coord[num_receptor_atoms][i];
}
/* Total up the partial charges as we go... */
q_tot += charge[num_receptor_atoms];
/* Increment the atom counter */
num_receptor_atoms++;
} /* endif */
} /* endwhile */
/* Finished reading in the lines of the receptor file */
(void) fclose( receptor_fileptr);
if ( has_receptor_types_in_gpf == 1 ) {
// Check that the number of atom types found in the receptor PDBQT
// file match the number parsed in by the "receptor_types" command
// in the GPF; if they do not match, exit!
if ( receptor_types_ct != receptor_types_gpf_ct ) {
(void) sprintf( message, "The number of atom types found in the receptor PDBQT (%d) does not match the number specified by the \"receptor_types\" command (%d) in the GPF!\n\n", receptor_types_ct, receptor_types_gpf_ct );
print_error( logFile, FATAL_ERROR, message );
// FATAL_ERROR will cause AutoGrid to exit...
}
}
/* Update the total number of atoms in the receptor */
(void) fprintf( logFile, "\nMaximum partial atomic charge found = %+.3lf e\n", q_max);