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