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