wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / initialize_p.F
1       block data
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.MCM'
5       include 'COMMON.MD'
6       data MovTypID
7      &  /'pool','chain regrow','multi-bond','phi','theta','side chain',
8      &   'total'/
9 c Conversion from poises to molecular unit and the gas constant
10       data cPoise /2.9361d0/, Rb /0.001986d0/
11       end
12 c--------------------------------------------------------------------------
13       subroutine initialize
14
15 C Define constants and zero out tables.
16 C
17       implicit real*8 (a-h,o-z)
18       include 'DIMENSIONS'
19 #ifdef MPI
20       include 'mpif.h'
21 #endif
22 #ifndef ISNAN
23       external proc_proc
24 #ifdef WINPGI
25 cMS$ATTRIBUTES C ::  proc_proc
26 #endif
27 #endif
28       include 'COMMON.IOUNITS'
29       include 'COMMON.CHAIN'
30       include 'COMMON.INTERACT'
31       include 'COMMON.GEO'
32       include 'COMMON.LOCAL'
33       include 'COMMON.TORSION'
34       include 'COMMON.FFIELD'
35       include 'COMMON.SBRIDGE'
36       include 'COMMON.MCM'
37       include 'COMMON.MINIM' 
38       include 'COMMON.DERIV'
39       include 'COMMON.SPLITELE'
40 c Common blocks from the diagonalization routines
41       COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
42       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
43       logical mask_r
44 c      real*8 text1 /'initial_i'/
45
46       mask_r=.false.
47 #ifndef ISNAN
48 c NaNQ initialization
49       i=-1
50       rr=dacos(100.0d0)
51 #ifdef WINPGI
52       idumm=proc_proc(rr,i)
53 #else
54       call proc_proc(rr,i)
55 #endif
56 #endif
57
58       kdiag=0
59       icorfl=0
60       iw=2
61 C
62 C The following is just to define auxiliary variables used in angle conversion
63 C
64       pi=4.0D0*datan(1.0D0)
65       dwapi=2.0D0*pi
66       dwapi3=dwapi/3.0D0
67       pipol=0.5D0*pi
68       deg2rad=pi/180.0D0
69       rad2deg=1.0D0/deg2rad
70       angmin=10.0D0*deg2rad
71 C
72 C Define I/O units.
73 C
74       inp=    1
75       iout=   2
76       ipdbin= 3
77       ipdb=   7
78       icart = 30
79       imol2=  4
80       igeom=  8
81       intin=  9
82       ithep= 11
83       ithep_pdb=51
84       irotam=12
85       irotam_pdb=52
86       itorp= 13
87       itordp= 23
88       ielep= 14
89       isidep=15 
90       iscpp=25
91       icbase=16
92       ifourier=20
93       istat= 17
94       irest1=55
95       irest2=56
96       iifrag=57
97       ientin=18
98       ientout=19
99       ibond = 28
100       isccor = 29
101 crc for write_rmsbank1  
102       izs1=21
103 cdr  include secondary structure prediction bias
104       isecpred=27
105 C
106 C CSA I/O units (separated from others especially for Jooyoung)
107 C
108       icsa_rbank=30
109       icsa_seed=31
110       icsa_history=32
111       icsa_bank=33
112       icsa_bank1=34
113       icsa_alpha=35
114       icsa_alpha1=36
115       icsa_bankt=37
116       icsa_int=39
117       icsa_bank_reminimized=38
118       icsa_native_int=41
119       icsa_in=40
120 crc for ifc error 118
121       icsa_pdb=42
122 C Lipidic input file for parameters range 60-79
123       iliptranpar=60
124 C input file for transfer sidechain and peptide group inside the
125 C lipidic environment if lipid is implicite
126
127 C DNA input files for parameters range 80-99
128 C Suger input files for parameters range 100-119
129 C All-atom input files for parameters range 120-149
130 C
131 C Set default weights of the energy terms.
132 C
133       wlong=1.0D0
134       welec=1.0D0
135       wtor =1.0D0
136       wang =1.0D0
137       wscloc=1.0D0
138       wstrain=1.0D0
139 C
140 C Zero out tables.
141 C
142 c      print '(a,$)','Inside initialize'
143 c      call memmon_print_usage()
144       do i=1,maxres2
145         do j=1,3
146           c(j,i)=0.0D0
147           dc(j,i)=0.0D0
148         enddo
149       enddo
150       do i=1,maxres
151         do j=1,3
152           xloc(j,i)=0.0D0
153         enddo
154       enddo
155       do i=1,ntyp
156         do j=1,ntyp
157           aa_aq(i,j)=0.0D0
158           bb_aq(i,j)=0.0D0
159           aa_lip(i,j)=0.0D0
160           bb_lip(i,j)=0.0D0
161           augm(i,j)=0.0D0
162           sigma(i,j)=0.0D0
163           r0(i,j)=0.0D0
164           chi(i,j)=0.0D0
165         enddo
166         do j=1,2
167           bad(i,j)=0.0D0
168         enddo
169         chip(i)=0.0D0
170         alp(i)=0.0D0
171         sigma0(i)=0.0D0
172         sigii(i)=0.0D0
173         rr0(i)=0.0D0
174         a0thet(i)=0.0D0
175         do j=1,2
176          do ichir1=-1,1
177           do ichir2=-1,1
178           athet(j,i,ichir1,ichir2)=0.0D0
179           bthet(j,i,ichir1,ichir2)=0.0D0
180           enddo
181          enddo
182         enddo
183         do j=0,3
184           polthet(j,i)=0.0D0
185         enddo
186         do j=1,3
187           gthet(j,i)=0.0D0
188         enddo
189         theta0(i)=0.0D0
190         sig0(i)=0.0D0
191         sigc0(i)=0.0D0
192         do j=1,maxlob
193           bsc(j,i)=0.0D0
194           do k=1,3
195             censc(k,j,i)=0.0D0
196           enddo
197           do k=1,3
198             do l=1,3
199               gaussc(l,k,j,i)=0.0D0
200             enddo
201           enddo
202           nlob(i)=0
203         enddo
204       enddo
205       nlob(ntyp1)=0
206       dsc(ntyp1)=0.0D0
207       do i=-maxtor,maxtor
208         itortyp(i)=0
209 cc      write (iout,*) "TU DOCHODZE",i,itortyp(i)
210        do iblock=1,2
211         do j=-maxtor,maxtor
212           do k=1,maxterm
213             v1(k,j,i,iblock)=0.0D0
214             v2(k,j,i,iblock)=0.0D0
215           enddo
216         enddo
217       enddo
218       enddo
219       do iblock=1,2
220        do i=-maxtor,maxtor
221         do j=-maxtor,maxtor
222          do k=-maxtor,maxtor
223           do l=1,maxtermd_1
224             v1c(1,l,i,j,k,iblock)=0.0D0
225             v1s(1,l,i,j,k,iblock)=0.0D0
226             v1c(2,l,i,j,k,iblock)=0.0D0
227             v1s(2,l,i,j,k,iblock)=0.0D0
228           enddo !l
229           do l=1,maxtermd_2
230            do m=1,maxtermd_2
231             v2c(m,l,i,j,k,iblock)=0.0D0
232             v2s(m,l,i,j,k,iblock)=0.0D0
233            enddo !m
234           enddo !l
235         enddo !k
236        enddo !j
237       enddo !i
238       enddo !iblock
239
240       do i=1,maxres
241         itype(i)=0
242         itel(i)=0
243       enddo
244 C Initialize the bridge arrays
245       ns=0
246       nss=0 
247       nhpb=0
248       do i=1,maxss
249         iss(i)=0
250       enddo
251       do i=1,maxdim
252         dhpb(i)=0.0D0
253       enddo
254       do i=1,maxres
255         ihpb(i)=0
256         jhpb(i)=0
257       enddo
258 C Initialize correlation arrays
259       do i=1,maxres
260        do k=1,2
261         b1(k,i)=0.0
262         b2(k,i)=0.0
263         b1tilde(k,i)=0.0
264 c        b2tilde(k,i)=0.0
265         do j=1,2
266 C        CC(j,k,i)=0.0
267 C        Ctilde(j,k,i)=0.0
268 C        DD(j,k,i)=0.0
269 C        Dtilde(j,k,i)=0.0
270         EE(j,k,i)=0.0
271         enddo
272        enddo
273       enddo
274       do i=-maxtor,maxtor
275        do k=1,2
276         do j=1,2
277         CC(j,k,i)=0.0
278         Ctilde(j,k,i)=0.0
279         DD(j,k,i)=0.0
280         Dtilde(j,k,i)=0.0
281         enddo
282       enddo
283       enddo
284 C
285 C Initialize timing.
286 C
287       call set_timers
288 C
289 C Initialize variables used in minimization.
290 C   
291 c     maxfun=5000
292 c     maxit=2000
293       maxfun=500
294       maxit=200
295       tolf=1.0D-2
296       rtolf=5.0D-4
297
298 C Initialize the variables responsible for the mode of gradient storage.
299 C
300       nfl=0
301       icg=1
302 C
303 C Initialize constants used to split the energy into long- and short-range
304 C components
305 C
306 C      r_cut=2.0d0
307 C      rlamb=0.3d0
308 #ifndef SPLITELE
309       nprint_ene=nprint_ene-1
310 #endif
311       return
312       end
313 c-------------------------------------------------------------------------
314       block data nazwy
315       implicit real*8 (a-h,o-z)
316       include 'DIMENSIONS'
317       include 'COMMON.NAMES'
318       include 'COMMON.FFIELD'
319       data restyp /
320      &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
321      & 'DSG','DGN','DSN','DTH',
322      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
323      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
324      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
325      &'AIB','ABU','D'/
326       data onelet /
327      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
328      &'a','y','w','v','l','i','f','m','c','x',
329      &'C','M','F','I','L','V','W','Y','A','G','T',
330      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
331       data potname /'LJ','LJK','BP','GB','GBV'/
332       data ename /
333      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
334      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
335      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB ","EVDWPP ",
336      &   "ESTR ","EVDW2_14 ","UCONST ", "EDIHC ","ESCCOR",
337      &   "Eliptran","Eafmforce","Ehomology","ESAXS"/
338       data wname /
339      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
340      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
341      &   "WSTRAIN","WVDWPP","WBOND","SCAL14","WDIHCSN","    ","WSCCOR",
342      &    "Wliptran","WAFM ","EHOMO","WSAXS"/
343       data nprint_ene /25/
344       data print_order/1,2,3,11,12,13,14,4,5,6,7,8,9,10,19,18,15,17,16,
345      & 21,24,22,23,20,25/
346       end 
347 c---------------------------------------------------------------------------
348       subroutine init_int_table
349       implicit real*8 (a-h,o-z)
350       include 'DIMENSIONS'
351 #ifdef MPI
352       include 'mpif.h'
353       integer blocklengths(15),displs(15)
354 #endif
355       include 'COMMON.CONTROL'
356       include 'COMMON.SETUP'
357       include 'COMMON.CHAIN'
358       include 'COMMON.INTERACT'
359       include 'COMMON.LOCAL'
360       include 'COMMON.SBRIDGE'
361       include 'COMMON.TORCNSTR'
362       include 'COMMON.IOUNITS'
363       include 'COMMON.DERIV'
364       include 'COMMON.CONTACTS'
365       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
366      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
367      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
368      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
369      & ielend_all(maxres,0:max_fg_procs-1),
370      & ntask_cont_from_all(0:max_fg_procs-1),
371      & itask_cont_from_all(0:max_fg_procs-1,0:max_fg_procs-1),
372      & ntask_cont_to_all(0:max_fg_procs-1),
373      & itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1)
374       integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP
375       logical scheck,lprint,flag
376 #ifdef MPI
377       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
378      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
379 C... Determine the numbers of start and end SC-SC interaction 
380 C... to deal with by current processor.
381       do i=0,nfgtasks-1
382         itask_cont_from(i)=fg_rank
383         itask_cont_to(i)=fg_rank
384       enddo
385       lprint=.false.
386       if (lprint)
387      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
388       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
389       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
390       if (lprint)
391      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
392      &  ' absolute rank',MyRank,
393      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
394      &  ' my_sc_inde',my_sc_inde
395       ind_sctint=0
396       iatsc_s=0
397       iatsc_e=0
398 #endif
399 c      lprint=.false.
400       do i=1,maxres
401         nint_gr(i)=0
402         nscp_gr(i)=0
403         do j=1,maxint_gr
404           istart(i,1)=0
405           iend(i,1)=0
406           ielstart(i)=0
407           ielend(i)=0
408           iscpstart(i,1)=0
409           iscpend(i,1)=0    
410         enddo
411       enddo
412       ind_scint=0
413       ind_scint_old=0
414 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
415 cd   &   (ihpb(i),jhpb(i),i=1,nss)
416       do i=nnt,nct-1
417         scheck=.false.
418         if (dyn_ss) goto 10
419         do ii=1,nss
420           if (ihpb(ii).eq.i+nres) then
421             scheck=.true.
422             jj=jhpb(ii)-nres
423             goto 10
424           endif
425         enddo
426    10   continue
427 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
428         if (scheck) then
429           if (jj.eq.i+1) then
430 #ifdef MPI
431 c            write (iout,*) 'jj=i+1'
432             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
433      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
434 #else
435             nint_gr(i)=1
436             istart(i,1)=i+2
437             iend(i,1)=nct
438 #endif
439           else if (jj.eq.nct) then
440 #ifdef MPI
441 c            write (iout,*) 'jj=nct'
442             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
443      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
444 #else
445             nint_gr(i)=1
446             istart(i,1)=i+1
447             iend(i,1)=nct-1
448 #endif
449           else
450 #ifdef MPI
451             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
452      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
453             ii=nint_gr(i)+1
454             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
455      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
456 #else
457             nint_gr(i)=2
458             istart(i,1)=i+1
459             iend(i,1)=jj-1
460             istart(i,2)=jj+1
461             iend(i,2)=nct
462 #endif
463           endif
464         else
465 #ifdef MPI
466           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
467      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
468 #else
469           nint_gr(i)=1
470           istart(i,1)=i+1
471           iend(i,1)=nct
472           ind_scint=ind_scint+nct-i
473 #endif
474         endif
475 #ifdef MPI
476         ind_scint_old=ind_scint
477 #endif
478       enddo
479    12 continue
480 #ifndef MPI
481       iatsc_s=nnt
482       iatsc_e=nct-1
483 #endif
484 #ifdef MPI
485       if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,
486      &   ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
487 #endif
488       if (lprint) then
489       write (iout,'(a)') 'Interaction array:'
490       do i=iatsc_s,iatsc_e
491         write (iout,'(i3,2(2x,2i3))') 
492      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
493       enddo
494       endif
495       ispp=4
496 #ifdef MPI
497 C Now partition the electrostatic-interaction array
498       npept=nct-nnt
499       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
500       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
501       if (lprint)
502      & write (*,*) 'Processor',fg_rank,' CG group',kolor,
503      &  ' absolute rank',MyRank,
504      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
505      &               ' my_ele_inde',my_ele_inde
506       iatel_s=0
507       iatel_e=0
508       ind_eleint=0
509       ind_eleint_old=0
510       do i=nnt,nct-3
511         ijunk=0
512         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
513      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
514       enddo ! i 
515    13 continue
516       if (iatel_s.eq.0) iatel_s=1
517       nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
518 c      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
519       call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
520 c      write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
521 c     & " my_ele_inde_vdw",my_ele_inde_vdw
522       ind_eleint_vdw=0
523       ind_eleint_vdw_old=0
524       iatel_s_vdw=0
525       iatel_e_vdw=0
526       do i=nnt,nct-3
527         ijunk=0
528         call int_partition(ind_eleint_vdw,my_ele_inds_vdw,
529      &    my_ele_inde_vdw,i,
530      &    iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),
531      &    ielend_vdw(i),*15)
532 c        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
533 c     &   " ielend_vdw",ielend_vdw(i)
534       enddo ! i 
535       if (iatel_s_vdw.eq.0) iatel_s_vdw=1
536    15 continue
537 #else
538       iatel_s=nnt
539       iatel_e=nct-5
540       do i=iatel_s,iatel_e
541         ielstart(i)=i+4
542         ielend(i)=nct-1
543       enddo
544       iatel_s_vdw=nnt
545       iatel_e_vdw=nct-3
546       do i=iatel_s_vdw,iatel_e_vdw
547         ielstart_vdw(i)=i+2
548         ielend_vdw(i)=nct-1
549       enddo
550 #endif
551       if (lprint) then
552         write (*,'(a)') 'Processor',fg_rank,' CG group',kolor,
553      &  ' absolute rank',MyRank
554         write (iout,*) 'Electrostatic interaction array:'
555         do i=iatel_s,iatel_e
556           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
557         enddo
558       endif ! lprint
559 c     iscp=3
560       iscp=2
561 C Partition the SC-p interaction array
562 #ifdef MPI
563       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
564       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
565       if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,
566      &  ' absolute rank',myrank,
567      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
568      &               ' my_scp_inde',my_scp_inde
569       iatscp_s=0
570       iatscp_e=0
571       ind_scpint=0
572       ind_scpint_old=0
573       do i=nnt,nct-1
574         if (i.lt.nnt+iscp) then
575 cd        write (iout,*) 'i.le.nnt+iscp'
576           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
577      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
578      &      iscpend(i,1),*14)
579         else if (i.gt.nct-iscp) then
580 cd        write (iout,*) 'i.gt.nct-iscp'
581           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
582      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
583      &      iscpend(i,1),*14)
584         else
585           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
586      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
587      &      iscpend(i,1),*14)
588           ii=nscp_gr(i)+1
589           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
590      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
591      &      iscpend(i,ii),*14)
592         endif
593       enddo ! i
594    14 continue
595 #else
596       iatscp_s=nnt
597       iatscp_e=nct-1
598       do i=nnt,nct-1
599         if (i.lt.nnt+iscp) then
600           nscp_gr(i)=1
601           iscpstart(i,1)=i+iscp
602           iscpend(i,1)=nct
603         elseif (i.gt.nct-iscp) then
604           nscp_gr(i)=1
605           iscpstart(i,1)=nnt
606           iscpend(i,1)=i-iscp
607         else
608           nscp_gr(i)=2
609           iscpstart(i,1)=nnt
610           iscpend(i,1)=i-iscp
611           iscpstart(i,2)=i+iscp
612           iscpend(i,2)=nct
613         endif 
614       enddo ! i
615 #endif
616       if (lprint) then
617         write (iout,'(a)') 'SC-p interaction array:'
618         do i=iatscp_s,iatscp_e
619           write (iout,'(i3,2(2x,2i3))') 
620      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
621         enddo
622       endif ! lprint
623 C Partition local interactions
624 #ifdef MPI
625       call int_bounds(nres-2,loc_start,loc_end)
626       loc_start=loc_start+1
627       loc_end=loc_end+1
628       call int_bounds(nres-2,ithet_start,ithet_end)
629       ithet_start=ithet_start+2
630       ithet_end=ithet_end+2
631       call int_bounds(nsaxs,isaxs_start,isaxs_end)
632 c      isaxs_start=isaxs_start+nnt-1
633 c      isaxs_end=isaxs_end+nnt-1
634       write (iout,*) me," isaxs_start",isaxs_start,
635      &  " isaxs_end",isaxs_end
636       call int_bounds(nct-nnt-2,iturn3_start,iturn3_end) 
637       iturn3_start=iturn3_start+nnt
638       iphi_start=iturn3_start+2
639       iturn3_end=iturn3_end+nnt
640       iphi_end=iturn3_end+2
641       iturn3_start=iturn3_start-1
642       iturn3_end=iturn3_end-1
643       call int_bounds(nres-3,itau_start,itau_end)
644       itau_start=itau_start+3
645       itau_end=itau_end+3
646       call int_bounds(nres-3,iphi1_start,iphi1_end)
647       iphi1_start=iphi1_start+3
648       iphi1_end=iphi1_end+3
649       call int_bounds(nct-nnt-3,iturn4_start,iturn4_end) 
650       iturn4_start=iturn4_start+nnt
651       iphid_start=iturn4_start+2
652       iturn4_end=iturn4_end+nnt
653       iphid_end=iturn4_end+2
654       iturn4_start=iturn4_start-1
655       iturn4_end=iturn4_end-1
656       call int_bounds(nres-2,ibond_start,ibond_end) 
657       ibond_start=ibond_start+1
658       ibond_end=ibond_end+1
659       call int_bounds(nct-nnt,ibondp_start,ibondp_end) 
660       ibondp_start=ibondp_start+nnt
661       ibondp_end=ibondp_end+nnt
662       call int_bounds(nres,ilip_start,ilip_end)
663       ilip_start=ilip_start
664       call int_bounds1(nres-1,ivec_start,ivec_end) 
665 c      print *,"Processor",myrank,fg_rank,fg_rank1,
666 c     &  " ivec_start",ivec_start," ivec_end",ivec_end
667       iset_start=loc_start+2
668       iset_end=loc_end+2
669       if (ndih_constr.eq.0) then
670         idihconstr_start=1
671         idihconstr_end=0
672       else
673         call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
674       endif
675 c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
676 c      nlen=nres-nnt+1
677       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
678       nlen=nres-nnt+1
679       call int_bounds(nsumgrad,ngrad_start,ngrad_end)
680       igrad_start=((2*nlen+1)
681      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
682       jgrad_start(igrad_start)=
683      &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
684      &    +igrad_start
685       jgrad_end(igrad_start)=nres
686       igrad_end=((2*nlen+1)
687      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
688       if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
689       jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
690      &    +igrad_end
691       do i=igrad_start+1,igrad_end-1
692         jgrad_start(i)=i+1
693         jgrad_end(i)=nres
694       enddo
695       if (lprint) then 
696         write (*,*) 'Processor:',fg_rank,' CG group',kolor,
697      & ' absolute rank',myrank,
698      & ' loc_start',loc_start,' loc_end',loc_end,
699      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
700      & ' iphi_start',iphi_start,' iphi_end',iphi_end,
701      & ' iphid_start',iphid_start,' iphid_end',iphid_end,
702      & ' ibond_start',ibond_start,' ibond_end',ibond_end,
703      & ' ibondp_start',ibondp_start,' ibondp_end',ibondp_end,
704      & ' iturn3_start',iturn3_start,' iturn3_end',iturn3_end,
705      & ' iturn4_start',iturn4_start,' iturn4_end',iturn4_end,
706      & ' ivec_start',ivec_start,' ivec_end',ivec_end,
707      & ' iset_start',iset_start,' iset_end',iset_end,
708      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
709      &   idihconstr_end
710        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
711      &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
712      &   ' ngrad_end',ngrad_end
713        do i=igrad_start,igrad_end
714          write(*,*) 'Processor:',fg_rank,myrank,i,
715      &    jgrad_start(i),jgrad_end(i)
716        enddo
717       endif
718       if (nfgtasks.gt.1) then
719         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
720      &    MPI_INTEGER,FG_COMM1,IERROR)
721         iaux=ivec_end-ivec_start+1
722         call MPI_Allgather(iaux,1,MPI_INTEGER,ivec_count(0),1,
723      &    MPI_INTEGER,FG_COMM1,IERROR)
724         call MPI_Allgather(iset_start-2,1,MPI_INTEGER,iset_displ(0),1,
725      &    MPI_INTEGER,FG_COMM,IERROR)
726         iaux=iset_end-iset_start+1
727         call MPI_Allgather(iaux,1,MPI_INTEGER,iset_count(0),1,
728      &    MPI_INTEGER,FG_COMM,IERROR)
729         call MPI_Allgather(ibond_start,1,MPI_INTEGER,ibond_displ(0),1,
730      &    MPI_INTEGER,FG_COMM,IERROR)
731         iaux=ibond_end-ibond_start+1
732         call MPI_Allgather(iaux,1,MPI_INTEGER,ibond_count(0),1,
733      &    MPI_INTEGER,FG_COMM,IERROR)
734         call MPI_Allgather(ithet_start,1,MPI_INTEGER,ithet_displ(0),1,
735      &    MPI_INTEGER,FG_COMM,IERROR)
736         iaux=ithet_end-ithet_start+1
737         call MPI_Allgather(iaux,1,MPI_INTEGER,ithet_count(0),1,
738      &    MPI_INTEGER,FG_COMM,IERROR)
739         call MPI_Allgather(iphi_start,1,MPI_INTEGER,iphi_displ(0),1,
740      &    MPI_INTEGER,FG_COMM,IERROR)
741         iaux=iphi_end-iphi_start+1
742         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi_count(0),1,
743      &    MPI_INTEGER,FG_COMM,IERROR)
744         call MPI_Allgather(iphi1_start,1,MPI_INTEGER,iphi1_displ(0),1,
745      &    MPI_INTEGER,FG_COMM,IERROR)
746         iaux=iphi1_end-iphi1_start+1
747         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi1_count(0),1,
748      &    MPI_INTEGER,FG_COMM,IERROR)
749         write (iout,*) "ielstart before zeroing out",max_fg_procs
750         call flush (iout)
751         do i=0,max_fg_procs-1
752           do j=1,maxres
753             ielstart_all(j,i)=0
754             ielend_all(j,i)=0
755           enddo
756         enddo
757         write (iout,*) "ielstart zeroed out"
758         call flush (iout)
759         call MPI_Allgather(iturn3_start,1,MPI_INTEGER,
760      &    iturn3_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
761         call MPI_Allgather(iturn4_start,1,MPI_INTEGER,
762      &    iturn4_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
763         call MPI_Allgather(iturn3_end,1,MPI_INTEGER,
764      &    iturn3_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
765         call MPI_Allgather(iturn4_end,1,MPI_INTEGER,
766      &    iturn4_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
767         call MPI_Allgather(iatel_s,1,MPI_INTEGER,
768      &    iatel_s_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
769         call MPI_Allgather(iatel_e,1,MPI_INTEGER,
770      &    iatel_e_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
771         call MPI_Allgather(ielstart(1),maxres,MPI_INTEGER,
772      &    ielstart_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
773         call MPI_Allgather(ielend(1),maxres,MPI_INTEGER,
774      &    ielend_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
775         if (lprint) then
776         write (iout,*) "iatel_s_all",(iatel_s_all(i),i=0,nfgtasks)
777         write (iout,*) "iatel_e_all",(iatel_e_all(i),i=0,nfgtasks)
778         write (iout,*) "iturn3_start_all",
779      &    (iturn3_start_all(i),i=0,nfgtasks-1)
780         write (iout,*) "iturn3_end_all",
781      &    (iturn3_end_all(i),i=0,nfgtasks-1)
782         write (iout,*) "iturn4_start_all",
783      &    (iturn4_start_all(i),i=0,nfgtasks-1)
784         write (iout,*) "iturn4_end_all",
785      &    (iturn4_end_all(i),i=0,nfgtasks-1)
786         write (iout,*) "The ielstart_all array"
787         do i=nnt,nct
788           write (iout,'(20i4)') i,(ielstart_all(i,j),j=0,nfgtasks-1)
789         enddo
790         write (iout,*) "The ielend_all array"
791         do i=nnt,nct
792           write (iout,'(20i4)') i,(ielend_all(i,j),j=0,nfgtasks-1)
793         enddo
794         call flush(iout)
795         endif
796         ntask_cont_from=0
797         ntask_cont_to=0
798         itask_cont_from(0)=fg_rank
799         itask_cont_to(0)=fg_rank
800         flag=.false.
801         do ii=iturn3_start,iturn3_end
802           call add_int(ii,ii+2,iturn3_sent(1,ii),
803      &                 ntask_cont_to,itask_cont_to,flag)
804         enddo
805         do ii=iturn4_start,iturn4_end
806           call add_int(ii,ii+3,iturn4_sent(1,ii),
807      &                 ntask_cont_to,itask_cont_to,flag)
808         enddo
809         do ii=iturn3_start,iturn3_end
810           call add_int_from(ii,ii+2,ntask_cont_from,itask_cont_from)
811         enddo
812         do ii=iturn4_start,iturn4_end
813           call add_int_from(ii,ii+3,ntask_cont_from,itask_cont_from)
814         enddo
815         if (lprint) then
816         write (iout,*) "After turn3 ntask_cont_from",ntask_cont_from,
817      &   " ntask_cont_to",ntask_cont_to
818         write (iout,*) "itask_cont_from",
819      &    (itask_cont_from(i),i=1,ntask_cont_from)
820         write (iout,*) "itask_cont_to",
821      &    (itask_cont_to(i),i=1,ntask_cont_to)
822         call flush(iout)
823         endif
824 c        write (iout,*) "Loop forward"
825 c        call flush(iout)
826         do i=iatel_s,iatel_e
827 c          write (iout,*) "from loop i=",i
828 c          call flush(iout)
829           do j=ielstart(i),ielend(i)
830             call add_int_from(i,j,ntask_cont_from,itask_cont_from)
831           enddo
832         enddo
833 c        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
834 c     &     " iatel_e",iatel_e
835 c        call flush(iout)
836         nat_sent=0
837         do i=iatel_s,iatel_e
838 c          write (iout,*) "i",i," ielstart",ielstart(i),
839 c     &      " ielend",ielend(i)
840 c          call flush(iout)
841           flag=.false.
842           do j=ielstart(i),ielend(i)
843             call add_int(i,j,iint_sent(1,j,nat_sent+1),ntask_cont_to,
844      &                   itask_cont_to,flag)
845           enddo
846           if (flag) then
847             nat_sent=nat_sent+1
848             iat_sent(nat_sent)=i
849           endif
850         enddo
851         if (lprint) then
852         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,
853      &   " ntask_cont_to",ntask_cont_to
854         write (iout,*) "itask_cont_from",
855      &    (itask_cont_from(i),i=1,ntask_cont_from)
856         write (iout,*) "itask_cont_to",
857      &    (itask_cont_to(i),i=1,ntask_cont_to)
858         call flush(iout)
859         write (iout,*) "iint_sent"
860         do i=1,nat_sent
861           ii=iat_sent(i)
862           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),
863      &      j=ielstart(ii),ielend(ii))
864         enddo
865         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,
866      &    " iturn3_end",iturn3_end
867         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),
868      &      i=iturn3_start,iturn3_end)
869         write (iout,*) "iturn4_sent iturn4_start",iturn4_start,
870      &    " iturn4_end",iturn4_end
871         write (iout,'(20i4)') (i,(iturn4_sent(j,i),j=1,4),
872      &      i=iturn4_start,iturn4_end)
873         call flush(iout)
874         endif
875         call MPI_Gather(ntask_cont_from,1,MPI_INTEGER,
876      &   ntask_cont_from_all,1,MPI_INTEGER,king,FG_COMM,IERR)
877 c        write (iout,*) "Gather ntask_cont_from ended"
878 c        call flush(iout)
879         call MPI_Gather(itask_cont_from(0),max_fg_procs,MPI_INTEGER,
880      &   itask_cont_from_all(0,0),max_fg_procs,MPI_INTEGER,king,
881      &   FG_COMM,IERR)
882 c        write (iout,*) "Gather itask_cont_from ended"
883 c        call flush(iout)
884         call MPI_Gather(ntask_cont_to,1,MPI_INTEGER,ntask_cont_to_all,
885      &   1,MPI_INTEGER,king,FG_COMM,IERR)
886 c        write (iout,*) "Gather ntask_cont_to ended"
887 c        call flush(iout)
888         call MPI_Gather(itask_cont_to,max_fg_procs,MPI_INTEGER,
889      &   itask_cont_to_all,max_fg_procs,MPI_INTEGER,king,FG_COMM,IERR)
890 c        write (iout,*) "Gather itask_cont_to ended"
891 c        call flush(iout)
892         if (fg_rank.eq.king) then
893          if (me.eq.0 .or. .not. out1file) then
894           write (iout,*)"Contact receive task map (proc, #tasks, tasks)"
895           do i=0,nfgtasks-1
896             write (iout,'(20i4)') i,ntask_cont_from_all(i),
897      &       (itask_cont_from_all(j,i),j=1,ntask_cont_from_all(i)) 
898           enddo
899           write (iout,*)
900           call flush(iout)
901           write (iout,*) "Contact send task map (proc, #tasks, tasks)"
902           do i=0,nfgtasks-1
903             write (iout,'(20i4)') i,ntask_cont_to_all(i),
904      &       (itask_cont_to_all(j,i),j=1,ntask_cont_to_all(i)) 
905           enddo
906           write (iout,*)
907           call flush(iout)
908          endif
909 C Check if every send will have a matching receive
910           ncheck_to=0
911           ncheck_from=0
912           do i=0,nfgtasks-1
913             ncheck_to=ncheck_to+ntask_cont_to_all(i)
914             ncheck_from=ncheck_from+ntask_cont_from_all(i)
915           enddo
916           if (me.eq.0 .or. .not. out1file)
917      &     write (iout,*) "Control sums",ncheck_from,ncheck_to
918           if (ncheck_from.ne.ncheck_to) then
919             write (iout,*) "Error: #receive differs from #send."
920             write (iout,*) "Terminating program...!"
921             call flush(iout)
922             flag=.false.
923           else
924             flag=.true.
925             do i=0,nfgtasks-1
926               do j=1,ntask_cont_to_all(i)
927                 ii=itask_cont_to_all(j,i)
928                 do k=1,ntask_cont_from_all(ii)
929                   if (itask_cont_from_all(k,ii).eq.i) then
930                     if(lprint)write(iout,*)"Matching send/receive",i,ii
931                     exit
932                   endif
933                 enddo
934                 if (k.eq.ntask_cont_from_all(ii)+1) then
935                   flag=.false.
936                   write (iout,*) "Error: send by",j," to",ii,
937      &            " would have no matching receive"
938                 endif
939               enddo
940             enddo
941           endif
942           if (.not.flag) then
943             write (iout,*) "Unmatched sends; terminating program"
944             call flush(iout)
945           endif
946         endif
947         call MPI_Bcast(flag,1,MPI_LOGICAL,king,FG_COMM,IERROR)
948 c        write (iout,*) "flag broadcast ended flag=",flag
949 c        call flush(iout)
950         if (.not.flag) then
951           call MPI_Finalize(IERROR)
952           stop "Error in INIT_INT_TABLE: unmatched send/receive."
953         endif
954         call MPI_Comm_group(FG_COMM,fg_group,IERR)
955 c        write (iout,*) "MPI_Comm_group ended"
956 c        call flush(iout)
957         call MPI_Group_incl(fg_group,ntask_cont_from+1,
958      &    itask_cont_from(0),CONT_FROM_GROUP,IERR)
959         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),
960      &    CONT_TO_GROUP,IERR)
961         do i=1,nat_sent
962           ii=iat_sent(i)
963           iaux=4*(ielend(ii)-ielstart(ii)+1)
964           call MPI_Group_translate_ranks(fg_group,iaux,
965      &      iint_sent(1,ielstart(ii),i),CONT_TO_GROUP, 
966      &      iint_sent_local(1,ielstart(ii),i),IERR )
967 c          write (iout,*) "Ranks translated i=",i
968 c          call flush(iout)
969         enddo
970         iaux=4*(iturn3_end-iturn3_start+1)
971         call MPI_Group_translate_ranks(fg_group,iaux,
972      &      iturn3_sent(1,iturn3_start),CONT_TO_GROUP,
973      &      iturn3_sent_local(1,iturn3_start),IERR)
974         iaux=4*(iturn4_end-iturn4_start+1)
975         call MPI_Group_translate_ranks(fg_group,iaux,
976      &      iturn4_sent(1,iturn4_start),CONT_TO_GROUP,
977      &      iturn4_sent_local(1,iturn4_start),IERR)
978         if (lprint) then
979         write (iout,*) "iint_sent_local"
980         do i=1,nat_sent
981           ii=iat_sent(i)
982           write (iout,'(20i4)') ii,(j,(iint_sent_local(k,j,i),k=1,4),
983      &      j=ielstart(ii),ielend(ii))
984           call flush(iout)
985         enddo
986         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,
987      &    " iturn3_end",iturn3_end
988         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),
989      &      i=iturn3_start,iturn3_end)
990         write (iout,*) "iturn4_sent_local iturn4_start",iturn4_start,
991      &    " iturn4_end",iturn4_end
992         write (iout,'(20i4)') (i,(iturn4_sent_local(j,i),j=1,4),
993      &      i=iturn4_start,iturn4_end)
994         call flush(iout)
995         endif
996         call MPI_Group_free(fg_group,ierr)
997         call MPI_Group_free(cont_from_group,ierr)
998         call MPI_Group_free(cont_to_group,ierr)
999         call MPI_Type_contiguous(3,MPI_DOUBLE_PRECISION,MPI_UYZ,IERROR)
1000         call MPI_Type_commit(MPI_UYZ,IERROR)
1001         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,
1002      &    IERROR)
1003         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
1004         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
1005         call MPI_Type_commit(MPI_MU,IERROR)
1006         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
1007         call MPI_Type_commit(MPI_MAT1,IERROR)
1008         call MPI_Type_contiguous(8,MPI_DOUBLE_PRECISION,MPI_MAT2,IERROR)
1009         call MPI_Type_commit(MPI_MAT2,IERROR)
1010         call MPI_Type_contiguous(6,MPI_DOUBLE_PRECISION,MPI_THET,IERROR)
1011         call MPI_Type_commit(MPI_THET,IERROR)
1012         call MPI_Type_contiguous(9,MPI_DOUBLE_PRECISION,MPI_GAM,IERROR)
1013         call MPI_Type_commit(MPI_GAM,IERROR)
1014 #ifndef MATGATHER
1015 c 9/22/08 Derived types to send matrices which appear in correlation terms
1016         do i=0,nfgtasks-1
1017           if (ivec_count(i).eq.ivec_count(0)) then
1018             lentyp(i)=0
1019           else
1020             lentyp(i)=1
1021           endif
1022         enddo
1023         do ind_typ=lentyp(0),lentyp(nfgtasks-1)
1024         if (ind_typ.eq.0) then
1025           ichunk=ivec_count(0)
1026         else
1027           ichunk=ivec_count(1)
1028         endif
1029 c        do i=1,4
1030 c          blocklengths(i)=4
1031 c        enddo
1032 c        displs(1)=0
1033 c        do i=2,4
1034 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1035 c        enddo
1036 c        do i=1,4
1037 c          blocklengths(i)=blocklengths(i)*ichunk
1038 c        enddo
1039 c        write (iout,*) "blocklengths and displs"
1040 c        do i=1,4
1041 c          write (iout,*) i,blocklengths(i),displs(i)
1042 c        enddo
1043 c        call flush(iout)
1044 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1045 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT1(ind_typ),IERROR)
1046 c        call MPI_Type_commit(MPI_ROTAT1(ind_typ),IERROR)
1047 c        write (iout,*) "MPI_ROTAT1",MPI_ROTAT1 
1048 c        do i=1,4
1049 c          blocklengths(i)=2
1050 c        enddo
1051 c        displs(1)=0
1052 c        do i=2,4
1053 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1054 c        enddo
1055 c        do i=1,4
1056 c          blocklengths(i)=blocklengths(i)*ichunk
1057 c        enddo
1058 c        write (iout,*) "blocklengths and displs"
1059 c        do i=1,4
1060 c          write (iout,*) i,blocklengths(i),displs(i)
1061 c        enddo
1062 c        call flush(iout)
1063 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1064 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT2(ind_typ),IERROR)
1065 c        call MPI_Type_commit(MPI_ROTAT2(ind_typ),IERROR)
1066 c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2 
1067         do i=1,8
1068           blocklengths(i)=2
1069         enddo
1070         displs(1)=0
1071         do i=2,8
1072           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1073         enddo
1074         do i=1,15
1075           blocklengths(i)=blocklengths(i)*ichunk
1076         enddo
1077         call MPI_Type_indexed(8,blocklengths,displs,
1078      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP11(ind_typ),IERROR)
1079         call MPI_Type_commit(MPI_PRECOMP11(ind_typ),IERROR)
1080         do i=1,8
1081           blocklengths(i)=4
1082         enddo
1083         displs(1)=0
1084         do i=2,8
1085           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1086         enddo
1087         do i=1,15
1088           blocklengths(i)=blocklengths(i)*ichunk
1089         enddo
1090         call MPI_Type_indexed(8,blocklengths,displs,
1091      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP12(ind_typ),IERROR)
1092         call MPI_Type_commit(MPI_PRECOMP12(ind_typ),IERROR)
1093         do i=1,6
1094           blocklengths(i)=4
1095         enddo
1096         displs(1)=0
1097         do i=2,6
1098           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1099         enddo
1100         do i=1,6
1101           blocklengths(i)=blocklengths(i)*ichunk
1102         enddo
1103         call MPI_Type_indexed(6,blocklengths,displs,
1104      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP22(ind_typ),IERROR)
1105         call MPI_Type_commit(MPI_PRECOMP22(ind_typ),IERROR)
1106         do i=1,2
1107           blocklengths(i)=8
1108         enddo
1109         displs(1)=0
1110         do i=2,2
1111           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1112         enddo
1113         do i=1,2
1114           blocklengths(i)=blocklengths(i)*ichunk
1115         enddo
1116         call MPI_Type_indexed(2,blocklengths,displs,
1117      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP23(ind_typ),IERROR)
1118         call MPI_Type_commit(MPI_PRECOMP23(ind_typ),IERROR)
1119         do i=1,4
1120           blocklengths(i)=1
1121         enddo
1122         displs(1)=0
1123         do i=2,4
1124           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1125         enddo
1126         do i=1,4
1127           blocklengths(i)=blocklengths(i)*ichunk
1128         enddo
1129         call MPI_Type_indexed(4,blocklengths,displs,
1130      &    MPI_DOUBLE_PRECISION,MPI_ROTAT_OLD(ind_typ),IERROR)
1131         call MPI_Type_commit(MPI_ROTAT_OLD(ind_typ),IERROR)
1132         enddo
1133 #endif
1134       endif
1135       iint_start=ivec_start+1
1136       iint_end=ivec_end+1
1137       do i=0,nfgtasks-1
1138           iint_count(i)=ivec_count(i)
1139           iint_displ(i)=ivec_displ(i)
1140           ivec_displ(i)=ivec_displ(i)-1
1141           iset_displ(i)=iset_displ(i)-1
1142           ithet_displ(i)=ithet_displ(i)-1
1143           iphi_displ(i)=iphi_displ(i)-1
1144           iphi1_displ(i)=iphi1_displ(i)-1
1145           ibond_displ(i)=ibond_displ(i)-1
1146       enddo
1147       if (nfgtasks.gt.1 .and. fg_rank.eq.king 
1148      &     .and. (me.eq.0 .or. .not. out1file)) then
1149         write (iout,*) "IVEC_DISPL, IVEC_COUNT, ISET_START, ISET_COUNT"
1150         do i=0,nfgtasks-1
1151           write (iout,*) i,ivec_displ(i),ivec_count(i),iset_displ(i),
1152      &      iset_count(i)
1153         enddo
1154         write (iout,*) "iphi_start",iphi_start," iphi_end",iphi_end,
1155      &    " iphi1_start",iphi1_start," iphi1_end",iphi1_end
1156         write (iout,*)"IPHI_COUNT, IPHI_DISPL, IPHI1_COUNT, IPHI1_DISPL"
1157         do i=0,nfgtasks-1
1158           write (iout,*) i,iphi_count(i),iphi_displ(i),iphi1_count(i),
1159      &      iphi1_displ(i)
1160         enddo
1161         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
1162      & nele_int_tot,' electrostatic and ',nscp_int_tot,
1163      & ' SC-p interactions','were distributed among',nfgtasks,
1164      & ' fine-grain processors.'
1165       endif
1166 #else
1167       loc_start=2
1168       loc_end=nres-1
1169       ithet_start=3 
1170       ithet_end=nres
1171       iturn3_start=nnt
1172       iturn3_end=nct-3
1173       iturn4_start=nnt
1174       iturn4_end=nct-4
1175       iphi_start=nnt+3
1176       iphi_end=nct
1177       iphi1_start=4
1178       iphi1_end=nres
1179       idihconstr_start=1
1180       idihconstr_end=ndih_constr
1181       isaxs_start=1
1182       isaxs_end=nsaxs
1183       iphid_start=iphi_start
1184       iphid_end=iphi_end-1
1185       itau_start=4
1186       itau_end=nres
1187       ibond_start=2
1188       ibond_end=nres-1
1189       ibondp_start=nnt
1190 C      ibondp_end=nct-1
1191       ibondp_end=nct
1192       ivec_start=1
1193       ivec_end=nres-1
1194       iset_start=3
1195       iset_end=nres+1
1196       iint_start=2
1197       iint_end=nres-1
1198       ilip_start=1
1199       ilip_end=nres
1200 #endif
1201       return
1202       end 
1203 #ifdef MPI
1204 c---------------------------------------------------------------------------
1205       subroutine add_int(ii,jj,itask,ntask_cont_to,itask_cont_to,flag)
1206       implicit none
1207       include "DIMENSIONS"
1208       include "COMMON.INTERACT"
1209       include "COMMON.SETUP"
1210       include "COMMON.IOUNITS"
1211       integer ii,jj,itask(4),ntask_cont_to,
1212      &itask_cont_to(0:max_fg_procs-1)
1213       logical flag
1214       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1215      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1216       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1217      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1218      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1219      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1220      & ielend_all(maxres,0:max_fg_procs-1)
1221       integer iproc,isent,k,l
1222 c Determines whether to send interaction ii,jj to other processors; a given
1223 c interaction can be sent to at most 2 processors.
1224 c Sets flag=.true. if interaction ii,jj needs to be sent to at least 
1225 c one processor, otherwise flag is unchanged from the input value.
1226       isent=0
1227       itask(1)=fg_rank
1228       itask(2)=fg_rank
1229       itask(3)=fg_rank
1230       itask(4)=fg_rank
1231 c      write (iout,*) "ii",ii," jj",jj
1232 c Loop over processors to check if anybody could need interaction ii,jj
1233       do iproc=0,fg_rank-1
1234 c Check if the interaction matches any turn3 at iproc
1235         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1236           l=k+2
1237           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1238      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1239      &    then 
1240 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj,"kl",k,l
1241 c            call flush(iout)
1242             flag=.true.
1243             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1244      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1245               isent=isent+1
1246               itask(isent)=iproc
1247               call add_task(iproc,ntask_cont_to,itask_cont_to)
1248             endif
1249           endif
1250         enddo
1251 C Check if the interaction matches any turn4 at iproc
1252         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1253           l=k+3
1254           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1255      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1256      &    then 
1257 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj," kl",k,l
1258 c            call flush(iout)
1259             flag=.true.
1260             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1261      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1262               isent=isent+1
1263               itask(isent)=iproc
1264               call add_task(iproc,ntask_cont_to,itask_cont_to)
1265             endif
1266           endif
1267         enddo
1268         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0 .and. 
1269      &  iatel_s_all(iproc).le.ii-1 .and. iatel_e_all(iproc).ge.ii-1)then
1270           if (ielstart_all(ii-1,iproc).le.jj-1.and.
1271      &        ielend_all(ii-1,iproc).ge.jj-1) then
1272             flag=.true.
1273             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1274      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1275               isent=isent+1
1276               itask(isent)=iproc
1277               call add_task(iproc,ntask_cont_to,itask_cont_to)
1278             endif
1279           endif
1280           if (ielstart_all(ii-1,iproc).le.jj+1.and.
1281      &        ielend_all(ii-1,iproc).ge.jj+1) then
1282             flag=.true.
1283             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1284      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1285               isent=isent+1
1286               itask(isent)=iproc
1287               call add_task(iproc,ntask_cont_to,itask_cont_to)
1288             endif
1289           endif
1290         endif
1291       enddo
1292       return
1293       end
1294 c---------------------------------------------------------------------------
1295       subroutine add_int_from(ii,jj,ntask_cont_from,itask_cont_from)
1296       implicit none
1297       include "DIMENSIONS"
1298       include "COMMON.INTERACT"
1299       include "COMMON.SETUP"
1300       include "COMMON.IOUNITS"
1301       integer ii,jj,itask(2),ntask_cont_from,
1302      & itask_cont_from(0:max_fg_procs-1)
1303       logical flag
1304       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1305      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1306       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1307      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1308      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1309      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1310      & ielend_all(maxres,0:max_fg_procs-1)
1311       integer iproc,k,l
1312       do iproc=fg_rank+1,nfgtasks-1
1313         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1314           l=k+2
1315           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1 
1316      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1) 
1317      &    then
1318 c            write (iout,*)"turn3 from iproc",iproc," ij",ii,jj," kl",k,l
1319             call add_task(iproc,ntask_cont_from,itask_cont_from)
1320           endif
1321         enddo 
1322         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1323           l=k+3
1324           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1 
1325      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1) 
1326      &    then
1327 c            write (iout,*)"turn4 from iproc",iproc," ij",ii,jj," kl",k,l
1328             call add_task(iproc,ntask_cont_from,itask_cont_from)
1329           endif
1330         enddo 
1331         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0) then
1332           if (ii+1.ge.iatel_s_all(iproc).and.ii+1.le.iatel_e_all(iproc))
1333      &    then
1334             if (jj+1.ge.ielstart_all(ii+1,iproc).and.
1335      &          jj+1.le.ielend_all(ii+1,iproc)) then
1336               call add_task(iproc,ntask_cont_from,itask_cont_from)
1337             endif            
1338             if (jj-1.ge.ielstart_all(ii+1,iproc).and.
1339      &          jj-1.le.ielend_all(ii+1,iproc)) then
1340               call add_task(iproc,ntask_cont_from,itask_cont_from)
1341             endif
1342           endif
1343           if (ii-1.ge.iatel_s_all(iproc).and.ii-1.le.iatel_e_all(iproc))
1344      &    then
1345             if (jj-1.ge.ielstart_all(ii-1,iproc).and.
1346      &          jj-1.le.ielend_all(ii-1,iproc)) then
1347               call add_task(iproc,ntask_cont_from,itask_cont_from)
1348             endif
1349             if (jj+1.ge.ielstart_all(ii-1,iproc).and.
1350      &          jj+1.le.ielend_all(ii-1,iproc)) then
1351                call add_task(iproc,ntask_cont_from,itask_cont_from)
1352             endif
1353           endif
1354         endif
1355       enddo
1356       return
1357       end
1358 c---------------------------------------------------------------------------
1359       subroutine add_task(iproc,ntask_cont,itask_cont)
1360       implicit none
1361       include "DIMENSIONS"
1362       integer iproc,ntask_cont,itask_cont(0:max_fg_procs-1)
1363       integer ii
1364       do ii=1,ntask_cont
1365         if (itask_cont(ii).eq.iproc) return
1366       enddo
1367       ntask_cont=ntask_cont+1
1368       itask_cont(ntask_cont)=iproc
1369       return
1370       end
1371 c---------------------------------------------------------------------------
1372       subroutine int_bounds(total_ints,lower_bound,upper_bound)
1373       implicit real*8 (a-h,o-z)
1374       include 'DIMENSIONS'
1375       include 'mpif.h'
1376       include 'COMMON.SETUP'
1377       integer total_ints,lower_bound,upper_bound
1378       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1379       nint=total_ints/nfgtasks
1380       do i=1,nfgtasks
1381         int4proc(i-1)=nint
1382       enddo
1383       nexcess=total_ints-nint*nfgtasks
1384       do i=1,nexcess
1385         int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1386       enddo
1387       lower_bound=0
1388       do i=0,fg_rank-1
1389         lower_bound=lower_bound+int4proc(i)
1390       enddo 
1391       upper_bound=lower_bound+int4proc(fg_rank)
1392       lower_bound=lower_bound+1
1393       return
1394       end
1395 c---------------------------------------------------------------------------
1396       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1397       implicit real*8 (a-h,o-z)
1398       include 'DIMENSIONS'
1399       include 'mpif.h'
1400       include 'COMMON.SETUP'
1401       integer total_ints,lower_bound,upper_bound
1402       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1403       nint=total_ints/nfgtasks1
1404       do i=1,nfgtasks1
1405         int4proc(i-1)=nint
1406       enddo
1407       nexcess=total_ints-nint*nfgtasks1
1408       do i=1,nexcess
1409         int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1410       enddo
1411       lower_bound=0
1412       do i=0,fg_rank1-1
1413         lower_bound=lower_bound+int4proc(i)
1414       enddo 
1415       upper_bound=lower_bound+int4proc(fg_rank1)
1416       lower_bound=lower_bound+1
1417       return
1418       end
1419 c---------------------------------------------------------------------------
1420       subroutine int_partition(int_index,lower_index,upper_index,atom,
1421      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
1422       implicit real*8 (a-h,o-z)
1423       include 'DIMENSIONS'
1424       include 'COMMON.IOUNITS'
1425       integer int_index,lower_index,upper_index,atom,at_start,at_end,
1426      & first_atom,last_atom,int_gr,jat_start,jat_end
1427       logical lprn
1428       lprn=.false.
1429       if (lprn) write (iout,*) 'int_index=',int_index
1430       int_index_old=int_index
1431       int_index=int_index+last_atom-first_atom+1
1432       if (lprn) 
1433      &   write (iout,*) 'int_index=',int_index,
1434      &               ' int_index_old',int_index_old,
1435      &               ' lower_index=',lower_index,
1436      &               ' upper_index=',upper_index,
1437      &               ' atom=',atom,' first_atom=',first_atom,
1438      &               ' last_atom=',last_atom
1439       if (int_index.ge.lower_index) then
1440         int_gr=int_gr+1
1441         if (at_start.eq.0) then
1442           at_start=atom
1443           jat_start=first_atom-1+lower_index-int_index_old
1444         else
1445           jat_start=first_atom
1446         endif
1447         if (lprn) write (iout,*) 'jat_start',jat_start
1448         if (int_index.ge.upper_index) then
1449           at_end=atom
1450           jat_end=first_atom-1+upper_index-int_index_old
1451           return1
1452         else
1453           jat_end=last_atom
1454         endif
1455         if (lprn) write (iout,*) 'jat_end',jat_end
1456       endif
1457       return
1458       end
1459 #endif
1460 c------------------------------------------------------------------------------
1461       subroutine hpb_partition
1462       implicit real*8 (a-h,o-z)
1463       include 'DIMENSIONS'
1464 #ifdef MPI
1465       include 'mpif.h'
1466 #endif
1467       include 'COMMON.SBRIDGE'
1468       include 'COMMON.IOUNITS'
1469       include 'COMMON.SETUP'
1470       include 'COMMON.CONTROL'
1471 c      write(2,*)"hpb_partition: nhpb=",nhpb
1472 #ifdef MPI
1473       call int_bounds(nhpb,link_start,link_end)
1474       if (.not. out1file) 
1475      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1476      &  ' absolute rank',MyRank,
1477      &  ' nhpb',nhpb,' link_start=',link_start,
1478      &  ' link_end',link_end
1479 #else
1480       link_start=1
1481       link_end=nhpb
1482 #endif
1483 c      write(2,*)"hpb_partition: link_start=",nhpb," link_end=",link_end
1484       return
1485       end
1486 c------------------------------------------------------------------------------
1487       subroutine homology_partition
1488       implicit real*8 (a-h,o-z)
1489       include 'DIMENSIONS'
1490 #ifdef MPI
1491       include 'mpif.h'
1492 #endif
1493       include 'COMMON.SBRIDGE'
1494       include 'COMMON.IOUNITS'
1495       include 'COMMON.SETUP'
1496       include 'COMMON.CONTROL'
1497       include 'COMMON.MD'
1498       include 'COMMON.INTERACT'
1499 cd      write(iout,*)"homology_partition: lim_odl=",lim_odl,
1500 cd     &   " lim_dih",lim_dih
1501 #ifdef MPI
1502       if (me.eq.king .or. .not. out1file) write (iout,*) "MPI"
1503       call int_bounds(lim_odl,link_start_homo,link_end_homo)
1504       call int_bounds(lim_dih,idihconstr_start_homo,
1505      &  idihconstr_end_homo)
1506       idihconstr_start_homo=idihconstr_start_homo+nnt-1+3
1507       idihconstr_end_homo=idihconstr_end_homo+nnt-1+3
1508       if (me.eq.king .or. .not. out1file) 
1509      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1510      &  ' absolute rank',MyRank,
1511      &  ' lim_odl',lim_odl,' link_start=',link_start_homo,
1512      &  ' link_end',link_end_homo,' lim_dih',lim_dih,
1513      &  ' idihconstr_start_homo',idihconstr_start_homo,
1514      &  ' idihconstr_end_homo',idihconstr_end_homo
1515 #else
1516       write (iout,*) "Not MPI"
1517       link_start_homo=1
1518       link_end_homo=lim_odl
1519       idihconstr_start_homo=nnt+3
1520       idihconstr_end_homo=lim_dih+nnt-1+3
1521       write (iout,*) 
1522      &  ' lim_odl',lim_odl,' link_start=',link_start_homo,
1523      &  ' link_end',link_end_homo,' lim_dih',lim_dih,
1524      &  ' idihconstr_start_homo',idihconstr_start_homo,
1525      &  ' idihconstr_end_homo',idihconstr_end_homo
1526 #endif
1527       return
1528       end