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