bd617ee787eeb901aef7bbd09367c3f697532ca9
[unres.git] / initialize_p.F
1       subroutine initialize
2
3 C Define constants and zero out tables.
4 C
5       implicit real*8 (a-h,o-z)
6       include 'DIMENSIONS'
7       include 'DIMENSIONS.ZSCOPT'
8 #ifdef MPI
9       include 'mpif.h'
10 #endif
11       include 'COMMON.IOUNITS'
12       include 'COMMON.CHAIN'
13       include 'COMMON.INTERACT'
14       include 'COMMON.GEO'
15       include 'COMMON.LOCAL'
16       include 'COMMON.TORSION'
17       include 'COMMON.FFIELD'
18       include 'COMMON.SBRIDGE'
19       include 'COMMON.MINIM' 
20       include 'COMMON.DERIV'
21       include "COMMON.WEIGHTS"
22       include "COMMON.NAMES"
23       include "COMMON.TIME1"
24 C
25 C The following is just to define auxiliary variables used in angle conversion
26 C
27       pi=4.0D0*datan(1.0D0)
28       dwapi=2.0D0*pi
29       dwapi3=dwapi/3.0D0
30       pipol=0.5D0*pi
31       deg2rad=pi/180.0D0
32       rad2deg=1.0D0/deg2rad
33       angmin=10.0D0*deg2rad
34 C
35 C Define I/O units.
36 C
37       inp=    1
38       iout=   2
39       ipdbin= 3
40       ipdb=   7
41       imol2=  4
42       igeom=  8
43       intin=  9
44       ithep= 11
45       irotam=12
46       itorp= 13
47       itordp= 23
48       ielep= 14
49       isidep=15 
50       isidep1=22
51       iscpp=25
52       icbase=16
53       ifourier=20
54       istat= 17
55       ientin=18
56       ientout=19
57       ibond=28
58       isccor=29
59 C
60 C WHAM files
61 C
62       ihist=30
63       iweight=31
64       izsc=32
65 C
66 C Set default weights of the energy terms.
67 C
68       wlong=1.0D0
69       welec=1.0D0
70       wtor =1.0D0
71       wang =1.0D0
72       wscloc=1.0D0
73       wstrain=1.0D0
74 C
75 C Zero out tables.
76 C
77       ndih_constr=0
78       do i=1,maxres2
79         do j=1,3
80           c(j,i)=0.0D0
81           dc(j,i)=0.0D0
82         enddo
83       enddo
84       do i=1,maxres
85         do j=1,3
86           xloc(j,i)=0.0D0
87         enddo
88       enddo
89       do i=1,ntyp
90         do j=1,ntyp
91           aa(i,j)=0.0D0
92           bb(i,j)=0.0D0
93           augm(i,j)=0.0D0
94           sigma(i,j)=0.0D0
95           r0(i,j)=0.0D0
96           chi(i,j)=0.0D0
97         enddo
98         do j=1,2
99           bad(i,j)=0.0D0
100         enddo
101         chip(i)=0.0D0
102         alp(i)=0.0D0
103         sigma0(i)=0.0D0
104         sigii(i)=0.0D0
105         rr0(i)=0.0D0
106         a0thet(i)=0.0D0
107         do j=1,2
108           athet(j,i)=0.0D0
109           bthet(j,i)=0.0D0
110         enddo
111         do j=0,3
112           polthet(j,i)=0.0D0
113         enddo
114         do j=1,3
115           gthet(j,i)=0.0D0
116         enddo
117         theta0(i)=0.0D0
118         sig0(i)=0.0D0
119         sigc0(i)=0.0D0
120         do j=1,maxlob
121           bsc(j,i)=0.0D0
122           do k=1,3
123             censc(k,j,i)=0.0D0
124           enddo
125           do k=1,3
126             do l=1,3
127               gaussc(l,k,j,i)=0.0D0
128             enddo
129           enddo
130           nlob(i)=0
131         enddo
132       enddo
133       nlob(ntyp1)=0
134       dsc(ntyp1)=0.0D0
135       do i=1,maxtor
136         itortyp(i)=0
137         do j=1,maxtor
138           do k=1,maxterm
139             v1(k,j,i)=0.0D0
140             v2(k,j,i)=0.0D0
141           enddo
142         enddo
143       enddo
144       do i=1,maxres
145         itype(i)=0
146         itel(i)=0
147       enddo
148 C Initialize the bridge arrays
149       ns=0
150       nss=0 
151       nhpb=0
152       do i=1,maxss
153         iss(i)=0
154       enddo
155       do i=1,maxdim
156         dhpb(i)=0.0D0
157       enddo
158       do i=1,maxres
159         ihpb(i)=0
160         jhpb(i)=0
161       enddo
162       do i=1,maxres
163         dyn_ss_mask(i)=.false.
164       enddo
165 C
166 C Initialize timing.
167 C
168       call set_timers
169 C
170 C Initialize variables used in minimization.
171 C   
172 c     maxfun=5000
173 c     maxit=2000
174       maxfun=500
175       maxit=200
176       tolf=1.0D-2
177       rtolf=5.0D-4
178
179 C Initialize the variables responsible for the mode of gradient storage.
180 C
181       nfl=0
182       icg=1
183       do i=1,14
184         do j=1,14
185           if (print_order(i).eq.j) then
186             iw(print_order(i))=j
187             goto 1121
188           endif
189         enddo
190 1121    continue
191       enddo
192       calc_grad=.false.
193 C Set timers and counters for the respective routines
194       t_func = 0.0d0
195       t_grad = 0.0d0
196       t_fhel = 0.0d0
197       t_fbet = 0.0d0
198       t_ghel = 0.0d0
199       t_gbet = 0.0d0
200       t_viol = 0.0d0
201       t_gviol = 0.0d0
202       n_func = 0
203       n_grad = 0
204       n_fhel = 0
205       n_fbet = 0
206       n_ghel = 0
207       n_gbet = 0
208       n_viol = 0
209       n_gviol = 0
210       n_map = 0
211 #ifndef SPLITELE
212       nprint_ene=nprint_ene-1
213 #endif
214       return
215       end
216 c-------------------------------------------------------------------------
217       block data nazwy
218       implicit real*8 (a-h,o-z)
219       include 'DIMENSIONS'
220       include 'DIMENSIONS.ZSCOPT'
221       include 'COMMON.NAMES'
222       include 'COMMON.WEIGHTS'
223       include 'COMMON.FFIELD'
224       data restyp /
225      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
226      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
227       data onelet /
228      &'C','M','F','I','L','V','W','Y','A','G','T',
229      &'S','Q','N','E','D','H','R','K','P','X'/
230       data potname /'LJ','LJK','BP','GB','GBV'/
231       data ename / 
232      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
233      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
234      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
235      &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T"/
236       data wname /
237      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
238      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
239      &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC"/
240       data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,
241      &    1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0,
242      &    0.0d0,0.0/
243       data nprint_ene /21/
244       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
245      &  16,15,17,20,21/
246       end 
247 c---------------------------------------------------------------------------
248       subroutine init_int_table
249       implicit real*8 (a-h,o-z)
250       include 'DIMENSIONS'
251       include 'DIMENSIONS.ZSCOPT'
252 #ifdef MPI
253       include 'mpif.h'
254 #endif
255 #ifdef MP
256       include 'COMMON.INFO'
257 #endif
258       include 'COMMON.CHAIN'
259       include 'COMMON.INTERACT'
260       include 'COMMON.LOCAL'
261       include 'COMMON.SBRIDGE'
262       include 'COMMON.IOUNITS'
263       logical scheck,lprint
264 #ifdef MPL
265       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
266      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
267 C... Determine the numbers of start and end SC-SC interaction 
268 C... to deal with by current processor.
269       lprint=.false.
270       if (lprint)
271      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
272       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
273       MyRank=MyID-(MyGroup-1)*fgProcs
274       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
275       if (lprint)
276      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
277      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
278      &  ' my_sc_inde',my_sc_inde
279       ind_sctint=0
280       iatsc_s=0
281       iatsc_e=0
282 #endif
283       lprint=.false.
284       do i=1,maxres
285         nint_gr(i)=0
286         nscp_gr(i)=0
287         do j=1,maxint_gr
288           istart(i,1)=0
289           iend(i,1)=0
290           ielstart(i)=0
291           ielend(i)=0
292           iscpstart(i,1)=0
293           iscpend(i,1)=0    
294         enddo
295       enddo
296       ind_scint=0
297       ind_scint_old=0
298 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
299 cd   &   (ihpb(i),jhpb(i),i=1,nss)
300       do i=nnt,nct-1
301         scheck=.false.
302         if (dyn_ss) go to 10
303         do ii=1,nss
304           if (ihpb(ii).eq.i+nres) then
305             scheck=.true.
306             jj=jhpb(ii)-nres
307             goto 10
308           endif
309         enddo
310    10   continue
311 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
312         if (scheck) then
313           if (jj.eq.i+1) then
314 #ifdef MPL
315             write (iout,*) 'jj=i+1'
316             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
317      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
318 #else
319             nint_gr(i)=1
320             istart(i,1)=i+2
321             iend(i,1)=nct
322 #endif
323           else if (jj.eq.nct) then
324 #ifdef MPL
325             write (iout,*) 'jj=nct'
326             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
327      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
328 #else
329             nint_gr(i)=1
330             istart(i,1)=i+1
331             iend(i,1)=nct-1
332 #endif
333           else
334 #ifdef MPL
335             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
336      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
337             ii=nint_gr(i)+1
338             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
339      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
340 #else
341             nint_gr(i)=2
342             istart(i,1)=i+1
343             iend(i,1)=jj-1
344             istart(i,2)=jj+1
345             iend(i,2)=nct
346 #endif
347           endif
348         else
349 #ifdef MPL
350           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
351      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
352 #else
353           nint_gr(i)=1
354           istart(i,1)=i+1
355           iend(i,1)=nct
356           ind_scint=int_scint+nct-i
357 #endif
358         endif
359 #ifdef MPL
360         ind_scint_old=ind_scint
361 #endif
362       enddo
363    12 continue
364 #ifndef MPL
365       iatsc_s=nnt
366       iatsc_e=nct-1
367 #endif
368 #ifdef MPL
369       if (lprint) then
370         write (iout,*) 'Processor',MyID,' Group',MyGroup
371         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
372       endif
373 #endif
374       if (lprint) then
375       write (iout,'(a)') 'Interaction array:'
376       do i=iatsc_s,iatsc_e
377         write (iout,'(i3,2(2x,2i3))') 
378      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
379       enddo
380       endif
381       ispp=2
382 #ifdef MPL
383 C Now partition the electrostatic-interaction array
384       npept=nct-nnt
385       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
386       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
387       if (lprint)
388      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
389      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
390      &               ' my_ele_inde',my_ele_inde
391       iatel_s=0
392       iatel_e=0
393       ind_eleint=0
394       ind_eleint_old=0
395       do i=nnt,nct-3
396         ijunk=0
397         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
398      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
399       enddo ! i 
400    13 continue
401 #else
402       iatel_s=nnt
403       iatel_e=nct-3
404       do i=iatel_s,iatel_e
405         ielstart(i)=i+2
406         ielend(i)=nct-1
407       enddo
408 #endif
409       if (lprint) then
410         write (iout,'(a)') 'Electrostatic interaction array:'
411         do i=iatel_s,iatel_e
412           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
413         enddo
414       endif ! lprint
415 c     iscp=3
416       iscp=2
417 C Partition the SC-p interaction array
418 #ifdef MPL
419       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
420       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
421       if (lprint)
422      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
423      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
424      &               ' my_scp_inde',my_scp_inde
425       iatscp_s=0
426       iatscp_e=0
427       ind_scpint=0
428       ind_scpint_old=0
429       do i=nnt,nct-1
430         if (i.lt.nnt+iscp) then
431 cd        write (iout,*) 'i.le.nnt+iscp'
432           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
433      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
434      &      iscpend(i,1),*14)
435         else if (i.gt.nct-iscp) then
436 cd        write (iout,*) 'i.gt.nct-iscp'
437           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
438      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
439      &      iscpend(i,1),*14)
440         else
441           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
442      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
443      &      iscpend(i,1),*14)
444           ii=nscp_gr(i)+1
445           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
446      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
447      &      iscpend(i,ii),*14)
448         endif
449       enddo ! i
450    14 continue
451 #else
452       iatscp_s=nnt
453       iatscp_e=nct-1
454       do i=nnt,nct-1
455         if (i.lt.nnt+iscp) then
456           nscp_gr(i)=1
457           iscpstart(i,1)=i+iscp
458           iscpend(i,1)=nct
459         elseif (i.gt.nct-iscp) then
460           nscp_gr(i)=1
461           iscpstart(i,1)=nnt
462           iscpend(i,1)=i-iscp
463         else
464           nscp_gr(i)=2
465           iscpstart(i,1)=nnt
466           iscpend(i,1)=i-iscp
467           iscpstart(i,2)=i+iscp
468           iscpend(i,2)=nct
469         endif 
470       enddo ! i
471 #endif
472               if (lprint) then
473         write (iout,'(a)') 'SC-p interaction array:'
474         do i=iatscp_s,iatscp_e
475           write (iout,'(i3,2(2x,2i3))') 
476      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
477         enddo
478       endif ! lprint
479 C Partition local interactions
480 #ifdef MPL
481       call int_bounds(nres-2,loc_start,loc_end)
482       loc_start=loc_start+1
483       loc_end=loc_end+1
484       call int_bounds(nres-2,ithet_start,ithet_end)
485       ithet_start=ithet_start+2
486       ithet_end=ithet_end+2
487       call int_bounds(nct-nnt-2,iphi_start,iphi_end) 
488       iphi_start=iphi_start+nnt+2
489       iphi_end=iphi_end+nnt+2
490       call int_bounds(nres-3,itau_start,itau_end)
491       itau_start=itau_start+3
492       itau_end=itau_end+3
493       if (lprint) then 
494         write (iout,*) 'Processor:',MyID,
495      & ' loc_start',loc_start,' loc_end',loc_end,
496      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
497      & ' iphi_start',iphi_start,' iphi_end',iphi_end
498         write (*,*) 'Processor:',MyID,
499      & ' loc_start',loc_start,' loc_end',loc_end,
500      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
501      & ' iphi_start',iphi_start,' iphi_end',iphi_end
502       endif
503       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
504         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
505      & nele_int_tot,' electrostatic and ',nscp_int_tot,
506      & ' SC-p interactions','were distributed among',fgprocs,
507      & ' fine-grain processors.'
508       endif
509 #else
510       loc_start=2
511       loc_end=nres-1
512       ithet_start=3 
513       ithet_end=nres
514       iphi_start=nnt+3
515       iphi_end=nct
516       itau_start=4
517       itau_end=nres
518 #endif
519       return
520       end 
521 c---------------------------------------------------------------------------
522       subroutine int_partition(int_index,lower_index,upper_index,atom,
523      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
524       implicit real*8 (a-h,o-z)
525       include 'DIMENSIONS'
526       include 'COMMON.IOUNITS'
527       integer int_index,lower_index,upper_index,atom,at_start,at_end,
528      & first_atom,last_atom,int_gr,jat_start,jat_end
529       logical lprn
530       lprn=.false.
531       if (lprn) write (iout,*) 'int_index=',int_index
532       int_index_old=int_index
533       int_index=int_index+last_atom-first_atom+1
534       if (lprn) 
535      &   write (iout,*) 'int_index=',int_index,
536      &               ' int_index_old',int_index_old,
537      &               ' lower_index=',lower_index,
538      &               ' upper_index=',upper_index,
539      &               ' atom=',atom,' first_atom=',first_atom,
540      &               ' last_atom=',last_atom
541       if (int_index.ge.lower_index) then
542         int_gr=int_gr+1
543         if (at_start.eq.0) then
544           at_start=atom
545           jat_start=first_atom-1+lower_index-int_index_old
546         else
547           jat_start=first_atom
548         endif
549         if (lprn) write (iout,*) 'jat_start',jat_start
550         if (int_index.ge.upper_index) then
551           at_end=atom
552           jat_end=first_atom-1+upper_index-int_index_old
553           return1
554         else
555           jat_end=last_atom
556         endif
557         if (lprn) write (iout,*) 'jat_end',jat_end
558       endif
559       return
560       end
561 c------------------------------------------------------------------------------
562       subroutine hpb_partition
563       implicit real*8 (a-h,o-z)
564       include 'DIMENSIONS'
565       include 'COMMON.SBRIDGE'
566       include 'COMMON.IOUNITS'
567 #ifdef MPL
568       include 'COMMON.INFO'
569       call int_bounds(nhpb,link_start,link_end)
570 #else
571       link_start=1
572       link_end=nhpb
573 #endif
574 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
575 cd   &  ' nhpb',nhpb,' link_start=',link_start,
576 cd   &  ' link_end',link_end
577       return
578       end