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