da0d3f9019b0ccefe24e576f091b6454a20c823b
[unres.git] / 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 ", "      ","ESCCOR"/
337       data wname /
338      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
339      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
340      &   "WSTRAIN","WVDWPP","WBOND","SCAL14","     ","    ","WSCCOR"/
341       data nprint_ene /20/
342       data print_order/1,2,3,11,12,13,14,4,5,6,7,8,9,10,19,18,15,17,16,
343      & 21,0/
344       end 
345 c---------------------------------------------------------------------------
346       subroutine init_int_table
347       implicit real*8 (a-h,o-z)
348       include 'DIMENSIONS'
349 #ifdef MPI
350       include 'mpif.h'
351       integer blocklengths(15),displs(15)
352 #endif
353       include 'COMMON.CONTROL'
354       include 'COMMON.SETUP'
355       include 'COMMON.CHAIN'
356       include 'COMMON.INTERACT'
357       include 'COMMON.LOCAL'
358       include 'COMMON.SBRIDGE'
359       include 'COMMON.TORCNSTR'
360       include 'COMMON.IOUNITS'
361       include 'COMMON.DERIV'
362       include 'COMMON.CONTACTS'
363       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
364      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
365      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
366      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
367      & ielend_all(maxres,0:max_fg_procs-1),
368      & ntask_cont_from_all(0:max_fg_procs-1),
369      & itask_cont_from_all(0:max_fg_procs-1,0:max_fg_procs-1),
370      & ntask_cont_to_all(0:max_fg_procs-1),
371      & itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1)
372       integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP
373       logical scheck,lprint,flag
374 #ifdef MPI
375       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
376      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
377 C... Determine the numbers of start and end SC-SC interaction 
378 C... to deal with by current processor.
379       do i=0,nfgtasks-1
380         itask_cont_from(i)=fg_rank
381         itask_cont_to(i)=fg_rank
382       enddo
383       lprint=energy_dec
384       if (lprint)
385      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
386       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
387       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
388       if (lprint)
389      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
390      &  ' absolute rank',MyRank,
391      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
392      &  ' my_sc_inde',my_sc_inde
393       ind_sctint=0
394       iatsc_s=0
395       iatsc_e=0
396 #endif
397 c      lprint=.false.
398       do i=1,maxres
399         nint_gr(i)=0
400         nscp_gr(i)=0
401         do j=1,maxint_gr
402           istart(i,1)=0
403           iend(i,1)=0
404           ielstart(i)=0
405           ielend(i)=0
406           iscpstart(i,1)=0
407           iscpend(i,1)=0    
408         enddo
409       enddo
410       ind_scint=0
411       ind_scint_old=0
412 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
413 cd   &   (ihpb(i),jhpb(i),i=1,nss)
414       do i=nnt,nct-1
415         scheck=.false.
416         if (dyn_ss) goto 10
417         do ii=1,nss
418           if (ihpb(ii).eq.i+nres) then
419             scheck=.true.
420             jj=jhpb(ii)-nres
421             goto 10
422           endif
423         enddo
424    10   continue
425 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
426         if (scheck) then
427           if (jj.eq.i+1) then
428 #ifdef MPI
429 c            write (iout,*) 'jj=i+1'
430             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
431      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
432 #else
433             nint_gr(i)=1
434             istart(i,1)=i+2
435             iend(i,1)=nct
436 #endif
437           else if (jj.eq.nct) then
438 #ifdef MPI
439 c            write (iout,*) 'jj=nct'
440             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
441      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
442 #else
443             nint_gr(i)=1
444             istart(i,1)=i+1
445             iend(i,1)=nct-1
446 #endif
447           else
448 #ifdef MPI
449             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
450      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
451             ii=nint_gr(i)+1
452             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
453      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
454 #else
455             nint_gr(i)=2
456             istart(i,1)=i+1
457             iend(i,1)=jj-1
458             istart(i,2)=jj+1
459             iend(i,2)=nct
460 #endif
461           endif
462         else
463 #ifdef MPI
464           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
465      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
466 #else
467           nint_gr(i)=1
468           istart(i,1)=i+1
469           iend(i,1)=nct
470           ind_scint=ind_scint+nct-i
471 #endif
472         endif
473 #ifdef MPI
474         ind_scint_old=ind_scint
475 #endif
476       enddo
477    12 continue
478 #ifndef MPI
479       iatsc_s=nnt
480       iatsc_e=nct-1
481 #endif
482       if (iatsc_s.eq.0) iatsc_s=1
483 #ifdef MPI
484       if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,
485      &   ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
486 #endif
487       if (lprint) then
488       write (iout,'(a)') 'Interaction array:'
489       do i=iatsc_s,iatsc_e
490         write (iout,'(i3,2(2x,2i3))') 
491      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
492       enddo
493       endif
494       ispp=4
495 #ifdef MPI
496 C Now partition the electrostatic-interaction array
497       npept=nct-nnt
498       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
499       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
500       if (lprint)
501      & write (*,*) 'Processor',fg_rank,' CG group',kolor,
502      &  ' absolute rank',MyRank,
503      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
504      &               ' my_ele_inde',my_ele_inde
505       iatel_s=0
506       iatel_e=0
507       ind_eleint=0
508       ind_eleint_old=0
509       do i=nnt,nct-3
510         ijunk=0
511         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
512      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
513       enddo ! i 
514    13 continue
515       if (iatel_s.eq.0) iatel_s=1
516       nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
517 c      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
518       call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
519 c      write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
520 c     & " my_ele_inde_vdw",my_ele_inde_vdw
521       ind_eleint_vdw=0
522       ind_eleint_vdw_old=0
523       iatel_s_vdw=0
524       iatel_e_vdw=0
525       do i=nnt,nct-3
526         ijunk=0
527         call int_partition(ind_eleint_vdw,my_ele_inds_vdw,
528      &    my_ele_inde_vdw,i,
529      &    iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),
530      &    ielend_vdw(i),*15)
531 c        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
532 c     &   " ielend_vdw",ielend_vdw(i)
533       enddo ! i 
534       if (iatel_s_vdw.eq.0) iatel_s_vdw=1
535    15 continue
536 #else
537       iatel_s=nnt
538       iatel_e=nct-5
539       do i=iatel_s,iatel_e
540         ielstart(i)=i+4
541         ielend(i)=nct-1
542       enddo
543       iatel_s_vdw=nnt
544       iatel_e_vdw=nct-3
545       do i=iatel_s_vdw,iatel_e_vdw
546         ielstart_vdw(i)=i+2
547         ielend_vdw(i)=nct-1
548       enddo
549 #endif
550       if (lprint) then
551         write (*,'(a)') 'Processor',fg_rank,' CG group',kolor,
552      &  ' absolute rank',MyRank
553         write (iout,*) 'Electrostatic interaction array:'
554         do i=iatel_s,iatel_e
555           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
556         enddo
557       endif ! lprint
558 c     iscp=3
559       iscp=2
560 C Partition the SC-p interaction array
561 #ifdef MPI
562       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
563       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
564       if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,
565      &  ' absolute rank',myrank,
566      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
567      &               ' my_scp_inde',my_scp_inde
568       iatscp_s=0
569       iatscp_e=0
570       ind_scpint=0
571       ind_scpint_old=0
572       do i=nnt,nct-1
573         if (i.lt.nnt+iscp) then
574 cd        write (iout,*) 'i.le.nnt+iscp'
575           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
576      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
577      &      iscpend(i,1),*14)
578         else if (i.gt.nct-iscp) then
579 cd        write (iout,*) 'i.gt.nct-iscp'
580           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
581      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
582      &      iscpend(i,1),*14)
583         else
584           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
585      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
586      &      iscpend(i,1),*14)
587           ii=nscp_gr(i)+1
588           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
589      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
590      &      iscpend(i,ii),*14)
591         endif
592       enddo ! i
593    14 continue
594 #else
595       iatscp_s=nnt
596       iatscp_e=nct-1
597       do i=nnt,nct-1
598         if (i.lt.nnt+iscp) then
599           nscp_gr(i)=1
600           iscpstart(i,1)=i+iscp
601           iscpend(i,1)=nct
602         elseif (i.gt.nct-iscp) then
603           nscp_gr(i)=1
604           iscpstart(i,1)=nnt
605           iscpend(i,1)=i-iscp
606         else
607           nscp_gr(i)=2
608           iscpstart(i,1)=nnt
609           iscpend(i,1)=i-iscp
610           iscpstart(i,2)=i+iscp
611           iscpend(i,2)=nct
612         endif 
613       enddo ! i
614 #endif
615       if (iatscp_s.eq.0) iatscp_s=1
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(nct-nnt-2,iturn3_start,iturn3_end) 
632       iturn3_start=iturn3_start+nnt
633       iphi_start=iturn3_start+2
634       iturn3_end=iturn3_end+nnt
635       iphi_end=iturn3_end+2
636       iturn3_start=iturn3_start-1
637       iturn3_end=iturn3_end-1
638       call int_bounds(nres-3,itau_start,itau_end)
639       itau_start=itau_start+3
640       itau_end=itau_end+3
641       call int_bounds(nres-3,iphi1_start,iphi1_end)
642       iphi1_start=iphi1_start+3
643       iphi1_end=iphi1_end+3
644       call int_bounds(nct-nnt-3,iturn4_start,iturn4_end) 
645       iturn4_start=iturn4_start+nnt
646       iphid_start=iturn4_start+2
647       iturn4_end=iturn4_end+nnt
648       iphid_end=iturn4_end+2
649       iturn4_start=iturn4_start-1
650       iturn4_end=iturn4_end-1
651       call int_bounds(nres-2,ibond_start,ibond_end) 
652       ibond_start=ibond_start+1
653       ibond_end=ibond_end+1
654       call int_bounds(nct-nnt,ibondp_start,ibondp_end) 
655       ibondp_start=ibondp_start+nnt
656       ibondp_end=ibondp_end+nnt
657       call int_bounds(nres,ilip_start,ilip_end)
658       ilip_start=ilip_start
659       call int_bounds1(nres-1,ivec_start,ivec_end) 
660 c      print *,"Processor",myrank,fg_rank,fg_rank1,
661 c     &  " ivec_start",ivec_start," ivec_end",ivec_end
662       iset_start=loc_start+2
663       iset_end=loc_end+2
664       if (ndih_constr.eq.0) then
665         idihconstr_start=1
666         idihconstr_end=0
667       else
668         call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
669       endif
670       if (ntheta_constr.eq.0) then
671         idihconstr_start=1
672         idihconstr_end=0
673       else
674         call int_bounds
675      &  (ntheta_constr,ithetaconstr_start,ithetaconstr_end)
676       endif
677 c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
678 c      nlen=nres-nnt+1
679       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
680       nlen=nres-nnt+1
681       call int_bounds(nsumgrad,ngrad_start,ngrad_end)
682       igrad_start=((2*nlen+1)
683      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
684       jgrad_start(igrad_start)=
685      &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
686      &    +igrad_start
687       jgrad_end(igrad_start)=nres
688       igrad_end=((2*nlen+1)
689      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
690       if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
691       jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
692      &    +igrad_end
693       do i=igrad_start+1,igrad_end-1
694         jgrad_start(i)=i+1
695         jgrad_end(i)=nres
696       enddo
697       if (lprint) then 
698         write (*,*) 'Processor:',fg_rank,' CG group',kolor,
699      & ' absolute rank',myrank,
700      & ' loc_start',loc_start,' loc_end',loc_end,
701      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
702      & ' iphi_start',iphi_start,' iphi_end',iphi_end,
703      & ' iphid_start',iphid_start,' iphid_end',iphid_end,
704      & ' ibond_start',ibond_start,' ibond_end',ibond_end,
705      & ' ibondp_start',ibondp_start,' ibondp_end',ibondp_end,
706      & ' iturn3_start',iturn3_start,' iturn3_end',iturn3_end,
707      & ' iturn4_start',iturn4_start,' iturn4_end',iturn4_end,
708      & ' ivec_start',ivec_start,' ivec_end',ivec_end,
709      & ' iset_start',iset_start,' iset_end',iset_end,
710      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
711      &   idihconstr_end,
712      & ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end',
713      &   ithetaconstr_end
714
715        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
716      &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
717      &   ' ngrad_end',ngrad_end
718        do i=igrad_start,igrad_end
719          write(*,*) 'Processor:',fg_rank,myrank,i,
720      &    jgrad_start(i),jgrad_end(i)
721        enddo
722       endif
723       if (nfgtasks.gt.1) then
724         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
725      &    MPI_INTEGER,FG_COMM1,IERROR)
726         iaux=ivec_end-ivec_start+1
727         call MPI_Allgather(iaux,1,MPI_INTEGER,ivec_count(0),1,
728      &    MPI_INTEGER,FG_COMM1,IERROR)
729         call MPI_Allgather(iset_start-2,1,MPI_INTEGER,iset_displ(0),1,
730      &    MPI_INTEGER,FG_COMM,IERROR)
731         iaux=iset_end-iset_start+1
732         call MPI_Allgather(iaux,1,MPI_INTEGER,iset_count(0),1,
733      &    MPI_INTEGER,FG_COMM,IERROR)
734         call MPI_Allgather(ibond_start,1,MPI_INTEGER,ibond_displ(0),1,
735      &    MPI_INTEGER,FG_COMM,IERROR)
736         iaux=ibond_end-ibond_start+1
737         call MPI_Allgather(iaux,1,MPI_INTEGER,ibond_count(0),1,
738      &    MPI_INTEGER,FG_COMM,IERROR)
739         call MPI_Allgather(ithet_start,1,MPI_INTEGER,ithet_displ(0),1,
740      &    MPI_INTEGER,FG_COMM,IERROR)
741         iaux=ithet_end-ithet_start+1
742         call MPI_Allgather(iaux,1,MPI_INTEGER,ithet_count(0),1,
743      &    MPI_INTEGER,FG_COMM,IERROR)
744         call MPI_Allgather(iphi_start,1,MPI_INTEGER,iphi_displ(0),1,
745      &    MPI_INTEGER,FG_COMM,IERROR)
746         iaux=iphi_end-iphi_start+1
747         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi_count(0),1,
748      &    MPI_INTEGER,FG_COMM,IERROR)
749         call MPI_Allgather(iphi1_start,1,MPI_INTEGER,iphi1_displ(0),1,
750      &    MPI_INTEGER,FG_COMM,IERROR)
751         iaux=iphi1_end-iphi1_start+1
752         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi1_count(0),1,
753      &    MPI_INTEGER,FG_COMM,IERROR)
754         do i=0,maxprocs-1
755           do j=1,maxres
756             ielstart_all(j,i)=0
757             ielend_all(j,i)=0
758           enddo
759         enddo
760         call MPI_Allgather(iturn3_start,1,MPI_INTEGER,
761      &    iturn3_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
762         call MPI_Allgather(iturn4_start,1,MPI_INTEGER,
763      &    iturn4_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
764         call MPI_Allgather(iturn3_end,1,MPI_INTEGER,
765      &    iturn3_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
766         call MPI_Allgather(iturn4_end,1,MPI_INTEGER,
767      &    iturn4_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
768         call MPI_Allgather(iatel_s,1,MPI_INTEGER,
769      &    iatel_s_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
770         call MPI_Allgather(iatel_e,1,MPI_INTEGER,
771      &    iatel_e_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
772         call MPI_Allgather(ielstart(1),maxres,MPI_INTEGER,
773      &    ielstart_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
774         call MPI_Allgather(ielend(1),maxres,MPI_INTEGER,
775      &    ielend_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
776         if (lprint) then
777         write (iout,*) "iatel_s_all",(iatel_s_all(i),i=0,nfgtasks)
778         write (iout,*) "iatel_e_all",(iatel_e_all(i),i=0,nfgtasks)
779         write (iout,*) "iturn3_start_all",
780      &    (iturn3_start_all(i),i=0,nfgtasks-1)
781         write (iout,*) "iturn3_end_all",
782      &    (iturn3_end_all(i),i=0,nfgtasks-1)
783         write (iout,*) "iturn4_start_all",
784      &    (iturn4_start_all(i),i=0,nfgtasks-1)
785         write (iout,*) "iturn4_end_all",
786      &    (iturn4_end_all(i),i=0,nfgtasks-1)
787         write (iout,*) "The ielstart_all array"
788         do i=nnt,nct
789           write (iout,'(20i4)') i,(ielstart_all(i,j),j=0,nfgtasks-1)
790         enddo
791         write (iout,*) "The ielend_all array"
792         do i=nnt,nct
793           write (iout,'(20i4)') i,(ielend_all(i,j),j=0,nfgtasks-1)
794         enddo
795         call flush(iout)
796         endif
797         ntask_cont_from=0
798         ntask_cont_to=0
799         itask_cont_from(0)=fg_rank
800         itask_cont_to(0)=fg_rank
801         flag=.false.
802         do ii=iturn3_start,iturn3_end
803           call add_int(ii,ii+2,iturn3_sent(1,ii),
804      &                 ntask_cont_to,itask_cont_to,flag)
805         enddo
806         do ii=iturn4_start,iturn4_end
807           call add_int(ii,ii+3,iturn4_sent(1,ii),
808      &                 ntask_cont_to,itask_cont_to,flag)
809         enddo
810         do ii=iturn3_start,iturn3_end
811           call add_int_from(ii,ii+2,ntask_cont_from,itask_cont_from)
812         enddo
813         do ii=iturn4_start,iturn4_end
814           call add_int_from(ii,ii+3,ntask_cont_from,itask_cont_from)
815         enddo
816         if (lprint) then
817         write (iout,*) "After turn3 ntask_cont_from",ntask_cont_from,
818      &   " ntask_cont_to",ntask_cont_to
819         write (iout,*) "itask_cont_from",
820      &    (itask_cont_from(i),i=1,ntask_cont_from)
821         write (iout,*) "itask_cont_to",
822      &    (itask_cont_to(i),i=1,ntask_cont_to)
823         call flush(iout)
824         endif
825 c        write (iout,*) "Loop forward"
826 c        call flush(iout)
827         do i=iatel_s,iatel_e
828 c          write (iout,*) "from loop i=",i
829 c          call flush(iout)
830           do j=ielstart(i),ielend(i)
831             call add_int_from(i,j,ntask_cont_from,itask_cont_from)
832           enddo
833         enddo
834 c        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
835 c     &     " iatel_e",iatel_e
836 c        call flush(iout)
837         nat_sent=0
838         do i=iatel_s,iatel_e
839 c          write (iout,*) "i",i," ielstart",ielstart(i),
840 c     &      " ielend",ielend(i)
841 c          call flush(iout)
842           flag=.false.
843           do j=ielstart(i),ielend(i)
844             call add_int(i,j,iint_sent(1,j,nat_sent+1),ntask_cont_to,
845      &                   itask_cont_to,flag)
846           enddo
847           if (flag) then
848             nat_sent=nat_sent+1
849             iat_sent(nat_sent)=i
850           endif
851         enddo
852         if (lprint) then
853         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,
854      &   " ntask_cont_to",ntask_cont_to
855         write (iout,*) "itask_cont_from",
856      &    (itask_cont_from(i),i=1,ntask_cont_from)
857         write (iout,*) "itask_cont_to",
858      &    (itask_cont_to(i),i=1,ntask_cont_to)
859         call flush(iout)
860         write (iout,*) "iint_sent"
861         do i=1,nat_sent
862           ii=iat_sent(i)
863           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),
864      &      j=ielstart(ii),ielend(ii))
865         enddo
866         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,
867      &    " iturn3_end",iturn3_end
868         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),
869      &      i=iturn3_start,iturn3_end)
870         write (iout,*) "iturn4_sent iturn4_start",iturn4_start,
871      &    " iturn4_end",iturn4_end
872         write (iout,'(20i4)') (i,(iturn4_sent(j,i),j=1,4),
873      &      i=iturn4_start,iturn4_end)
874         call flush(iout)
875         endif
876         call MPI_Gather(ntask_cont_from,1,MPI_INTEGER,
877      &   ntask_cont_from_all,1,MPI_INTEGER,king,FG_COMM,IERR)
878 c        write (iout,*) "Gather ntask_cont_from ended"
879 c        call flush(iout)
880         call MPI_Gather(itask_cont_from(0),max_fg_procs,MPI_INTEGER,
881      &   itask_cont_from_all(0,0),max_fg_procs,MPI_INTEGER,king,
882      &   FG_COMM,IERR)
883 c        write (iout,*) "Gather itask_cont_from ended"
884 c        call flush(iout)
885         call MPI_Gather(ntask_cont_to,1,MPI_INTEGER,ntask_cont_to_all,
886      &   1,MPI_INTEGER,king,FG_COMM,IERR)
887 c        write (iout,*) "Gather ntask_cont_to ended"
888 c        call flush(iout)
889         call MPI_Gather(itask_cont_to,max_fg_procs,MPI_INTEGER,
890      &   itask_cont_to_all,max_fg_procs,MPI_INTEGER,king,FG_COMM,IERR)
891 c        write (iout,*) "Gather itask_cont_to ended"
892 c        call flush(iout)
893         if (fg_rank.eq.king) 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 C Check if every send will have a matching receive
909           ncheck_to=0
910           ncheck_from=0
911           do i=0,nfgtasks-1
912             ncheck_to=ncheck_to+ntask_cont_to_all(i)
913             ncheck_from=ncheck_from+ntask_cont_from_all(i)
914           enddo
915           write (iout,*) "Control sums",ncheck_from,ncheck_to
916           if (ncheck_from.ne.ncheck_to) then
917             write (iout,*) "Error: #receive differs from #send."
918             write (iout,*) "Terminating program...!"
919             call flush(iout)
920             flag=.false.
921           else
922             flag=.true.
923             do i=0,nfgtasks-1
924               do j=1,ntask_cont_to_all(i)
925                 ii=itask_cont_to_all(j,i)
926                 do k=1,ntask_cont_from_all(ii)
927                   if (itask_cont_from_all(k,ii).eq.i) then
928                     if(lprint)write(iout,*)"Matching send/receive",i,ii
929                     exit
930                   endif
931                 enddo
932                 if (k.eq.ntask_cont_from_all(ii)+1) then
933                   flag=.false.
934                   write (iout,*) "Error: send by",j," to",ii,
935      &            " would have no matching receive"
936                 endif
937               enddo
938             enddo
939           endif
940           if (.not.flag) then
941             write (iout,*) "Unmatched sends; terminating program"
942             call flush(iout)
943           endif
944         endif
945         call MPI_Bcast(flag,1,MPI_LOGICAL,king,FG_COMM,IERROR)
946 c        write (iout,*) "flag broadcast ended flag=",flag
947 c        call flush(iout)
948         if (.not.flag) then
949           call MPI_Finalize(IERROR)
950           stop "Error in INIT_INT_TABLE: unmatched send/receive."
951         endif
952         call MPI_Comm_group(FG_COMM,fg_group,IERR)
953 c        write (iout,*) "MPI_Comm_group ended"
954 c        call flush(iout)
955         call MPI_Group_incl(fg_group,ntask_cont_from+1,
956      &    itask_cont_from(0),CONT_FROM_GROUP,IERR)
957         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),
958      &    CONT_TO_GROUP,IERR)
959         do i=1,nat_sent
960           ii=iat_sent(i)
961           iaux=4*(ielend(ii)-ielstart(ii)+1)
962           call MPI_Group_translate_ranks(fg_group,iaux,
963      &      iint_sent(1,ielstart(ii),i),CONT_TO_GROUP, 
964      &      iint_sent_local(1,ielstart(ii),i),IERR )
965 c          write (iout,*) "Ranks translated i=",i
966 c          call flush(iout)
967         enddo
968         iaux=4*(iturn3_end-iturn3_start+1)
969         call MPI_Group_translate_ranks(fg_group,iaux,
970      &      iturn3_sent(1,iturn3_start),CONT_TO_GROUP,
971      &      iturn3_sent_local(1,iturn3_start),IERR)
972         iaux=4*(iturn4_end-iturn4_start+1)
973         call MPI_Group_translate_ranks(fg_group,iaux,
974      &      iturn4_sent(1,iturn4_start),CONT_TO_GROUP,
975      &      iturn4_sent_local(1,iturn4_start),IERR)
976         if (lprint) then
977         write (iout,*) "iint_sent_local"
978         do i=1,nat_sent
979           ii=iat_sent(i)
980           write (iout,'(20i4)') ii,(j,(iint_sent_local(k,j,i),k=1,4),
981      &      j=ielstart(ii),ielend(ii))
982           call flush(iout)
983         enddo
984         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,
985      &    " iturn3_end",iturn3_end
986         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),
987      &      i=iturn3_start,iturn3_end)
988         write (iout,*) "iturn4_sent_local iturn4_start",iturn4_start,
989      &    " iturn4_end",iturn4_end
990         write (iout,'(20i4)') (i,(iturn4_sent_local(j,i),j=1,4),
991      &      i=iturn4_start,iturn4_end)
992         call flush(iout)
993         endif
994         call MPI_Group_free(fg_group,ierr)
995         call MPI_Group_free(cont_from_group,ierr)
996         call MPI_Group_free(cont_to_group,ierr)
997         call MPI_Type_contiguous(3,MPI_DOUBLE_PRECISION,MPI_UYZ,IERROR)
998         call MPI_Type_commit(MPI_UYZ,IERROR)
999         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,
1000      &    IERROR)
1001         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
1002         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
1003         call MPI_Type_commit(MPI_MU,IERROR)
1004         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
1005         call MPI_Type_commit(MPI_MAT1,IERROR)
1006         call MPI_Type_contiguous(8,MPI_DOUBLE_PRECISION,MPI_MAT2,IERROR)
1007         call MPI_Type_commit(MPI_MAT2,IERROR)
1008         call MPI_Type_contiguous(6,MPI_DOUBLE_PRECISION,MPI_THET,IERROR)
1009         call MPI_Type_commit(MPI_THET,IERROR)
1010         call MPI_Type_contiguous(9,MPI_DOUBLE_PRECISION,MPI_GAM,IERROR)
1011         call MPI_Type_commit(MPI_GAM,IERROR)
1012 #ifndef MATGATHER
1013 c 9/22/08 Derived types to send matrices which appear in correlation terms
1014         do i=0,nfgtasks-1
1015           if (ivec_count(i).eq.ivec_count(0)) then
1016             lentyp(i)=0
1017           else
1018             lentyp(i)=1
1019           endif
1020         enddo
1021         do ind_typ=lentyp(0),lentyp(nfgtasks-1)
1022         if (ind_typ.eq.0) then
1023           ichunk=ivec_count(0)
1024         else
1025           ichunk=ivec_count(1)
1026         endif
1027 c        do i=1,4
1028 c          blocklengths(i)=4
1029 c        enddo
1030 c        displs(1)=0
1031 c        do i=2,4
1032 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1033 c        enddo
1034 c        do i=1,4
1035 c          blocklengths(i)=blocklengths(i)*ichunk
1036 c        enddo
1037 c        write (iout,*) "blocklengths and displs"
1038 c        do i=1,4
1039 c          write (iout,*) i,blocklengths(i),displs(i)
1040 c        enddo
1041 c        call flush(iout)
1042 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1043 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT1(ind_typ),IERROR)
1044 c        call MPI_Type_commit(MPI_ROTAT1(ind_typ),IERROR)
1045 c        write (iout,*) "MPI_ROTAT1",MPI_ROTAT1 
1046 c        do i=1,4
1047 c          blocklengths(i)=2
1048 c        enddo
1049 c        displs(1)=0
1050 c        do i=2,4
1051 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1052 c        enddo
1053 c        do i=1,4
1054 c          blocklengths(i)=blocklengths(i)*ichunk
1055 c        enddo
1056 c        write (iout,*) "blocklengths and displs"
1057 c        do i=1,4
1058 c          write (iout,*) i,blocklengths(i),displs(i)
1059 c        enddo
1060 c        call flush(iout)
1061 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1062 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT2(ind_typ),IERROR)
1063 c        call MPI_Type_commit(MPI_ROTAT2(ind_typ),IERROR)
1064 c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2 
1065         do i=1,8
1066           blocklengths(i)=2
1067         enddo
1068         displs(1)=0
1069         do i=2,8
1070           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1071         enddo
1072         do i=1,15
1073           blocklengths(i)=blocklengths(i)*ichunk
1074         enddo
1075         call MPI_Type_indexed(8,blocklengths,displs,
1076      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP11(ind_typ),IERROR)
1077         call MPI_Type_commit(MPI_PRECOMP11(ind_typ),IERROR)
1078         do i=1,8
1079           blocklengths(i)=4
1080         enddo
1081         displs(1)=0
1082         do i=2,8
1083           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1084         enddo
1085         do i=1,15
1086           blocklengths(i)=blocklengths(i)*ichunk
1087         enddo
1088         call MPI_Type_indexed(8,blocklengths,displs,
1089      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP12(ind_typ),IERROR)
1090         call MPI_Type_commit(MPI_PRECOMP12(ind_typ),IERROR)
1091         do i=1,6
1092           blocklengths(i)=4
1093         enddo
1094         displs(1)=0
1095         do i=2,6
1096           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1097         enddo
1098         do i=1,6
1099           blocklengths(i)=blocklengths(i)*ichunk
1100         enddo
1101         call MPI_Type_indexed(6,blocklengths,displs,
1102      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP22(ind_typ),IERROR)
1103         call MPI_Type_commit(MPI_PRECOMP22(ind_typ),IERROR)
1104         do i=1,2
1105           blocklengths(i)=8
1106         enddo
1107         displs(1)=0
1108         do i=2,2
1109           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1110         enddo
1111         do i=1,2
1112           blocklengths(i)=blocklengths(i)*ichunk
1113         enddo
1114         call MPI_Type_indexed(2,blocklengths,displs,
1115      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP23(ind_typ),IERROR)
1116         call MPI_Type_commit(MPI_PRECOMP23(ind_typ),IERROR)
1117         do i=1,4
1118           blocklengths(i)=1
1119         enddo
1120         displs(1)=0
1121         do i=2,4
1122           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1123         enddo
1124         do i=1,4
1125           blocklengths(i)=blocklengths(i)*ichunk
1126         enddo
1127         call MPI_Type_indexed(4,blocklengths,displs,
1128      &    MPI_DOUBLE_PRECISION,MPI_ROTAT_OLD(ind_typ),IERROR)
1129         call MPI_Type_commit(MPI_ROTAT_OLD(ind_typ),IERROR)
1130         enddo
1131 #endif
1132       endif
1133       iint_start=ivec_start+1
1134       iint_end=ivec_end+1
1135       do i=0,nfgtasks-1
1136           iint_count(i)=ivec_count(i)
1137           iint_displ(i)=ivec_displ(i)
1138           ivec_displ(i)=ivec_displ(i)-1
1139           iset_displ(i)=iset_displ(i)-1
1140           ithet_displ(i)=ithet_displ(i)-1
1141           iphi_displ(i)=iphi_displ(i)-1
1142           iphi1_displ(i)=iphi1_displ(i)-1
1143           ibond_displ(i)=ibond_displ(i)-1
1144       enddo
1145       if (nfgtasks.gt.1 .and. fg_rank.eq.king 
1146      &     .and. (me.eq.0 .or. .not. out1file)) then
1147         write (iout,*) "IVEC_DISPL, IVEC_COUNT, ISET_START, ISET_COUNT"
1148         do i=0,nfgtasks-1
1149           write (iout,*) i,ivec_displ(i),ivec_count(i),iset_displ(i),
1150      &      iset_count(i)
1151         enddo
1152         write (iout,*) "iphi_start",iphi_start," iphi_end",iphi_end,
1153      &    " iphi1_start",iphi1_start," iphi1_end",iphi1_end
1154         write (iout,*)"IPHI_COUNT, IPHI_DISPL, IPHI1_COUNT, IPHI1_DISPL"
1155         do i=0,nfgtasks-1
1156           write (iout,*) i,iphi_count(i),iphi_displ(i),iphi1_count(i),
1157      &      iphi1_displ(i)
1158         enddo
1159         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
1160      & nele_int_tot,' electrostatic and ',nscp_int_tot,
1161      & ' SC-p interactions','were distributed among',nfgtasks,
1162      & ' fine-grain processors.'
1163       endif
1164 #else
1165       loc_start=2
1166       loc_end=nres-1
1167       ithet_start=3 
1168       ithet_end=nres
1169       iturn3_start=nnt
1170       iturn3_end=nct-3
1171       iturn4_start=nnt
1172       iturn4_end=nct-4
1173       iphi_start=nnt+3
1174       iphi_end=nct
1175       iphi1_start=4
1176       iphi1_end=nres
1177       idihconstr_start=1
1178       idihconstr_end=ndih_constr
1179       ithetaconstr_start=1
1180       ithetaconstr_end=ntheta_constr
1181       iphid_start=iphi_start
1182       iphid_end=iphi_end-1
1183       itau_start=4
1184       itau_end=nres
1185       ibond_start=2
1186       ibond_end=nres-1
1187       ibondp_start=nnt
1188 C      ibondp_end=nct-1
1189       ibondp_end=nct
1190       ivec_start=1
1191       ivec_end=nres-1
1192       iset_start=3
1193       iset_end=nres+1
1194       iint_start=2
1195       iint_end=nres-1
1196       ilip_start=1
1197       ilip_end=nres
1198 #endif
1199       return
1200       end 
1201 #ifdef MPI
1202 c---------------------------------------------------------------------------
1203       subroutine add_int(ii,jj,itask,ntask_cont_to,itask_cont_to,flag)
1204       implicit none
1205       include "DIMENSIONS"
1206       include "COMMON.INTERACT"
1207       include "COMMON.SETUP"
1208       include "COMMON.IOUNITS"
1209       integer ii,jj,itask(4),ntask_cont_to,
1210      &itask_cont_to(0:max_fg_procs-1)
1211       logical flag
1212       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1213      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1214       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1215      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1216      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1217      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1218      & ielend_all(maxres,0:max_fg_procs-1)
1219       integer iproc,isent,k,l
1220 c Determines whether to send interaction ii,jj to other processors; a given
1221 c interaction can be sent to at most 2 processors.
1222 c Sets flag=.true. if interaction ii,jj needs to be sent to at least 
1223 c one processor, otherwise flag is unchanged from the input value.
1224       isent=0
1225       itask(1)=fg_rank
1226       itask(2)=fg_rank
1227       itask(3)=fg_rank
1228       itask(4)=fg_rank
1229 c      write (iout,*) "ii",ii," jj",jj
1230 c Loop over processors to check if anybody could need interaction ii,jj
1231       do iproc=0,fg_rank-1
1232 c Check if the interaction matches any turn3 at iproc
1233         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1234           l=k+2
1235           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1236      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1237      &    then 
1238 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj,"kl",k,l
1239 c            call flush(iout)
1240             flag=.true.
1241             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1242      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1243               isent=isent+1
1244               itask(isent)=iproc
1245               call add_task(iproc,ntask_cont_to,itask_cont_to)
1246             endif
1247           endif
1248         enddo
1249 C Check if the interaction matches any turn4 at iproc
1250         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1251           l=k+3
1252           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1253      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1254      &    then 
1255 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj," kl",k,l
1256 c            call flush(iout)
1257             flag=.true.
1258             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1259      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1260               isent=isent+1
1261               itask(isent)=iproc
1262               call add_task(iproc,ntask_cont_to,itask_cont_to)
1263             endif
1264           endif
1265         enddo
1266         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0 .and. 
1267      &  iatel_s_all(iproc).le.ii-1 .and. iatel_e_all(iproc).ge.ii-1)then
1268           if (ielstart_all(ii-1,iproc).le.jj-1.and.
1269      &        ielend_all(ii-1,iproc).ge.jj-1) then
1270             flag=.true.
1271             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1272      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1273               isent=isent+1
1274               itask(isent)=iproc
1275               call add_task(iproc,ntask_cont_to,itask_cont_to)
1276             endif
1277           endif
1278           if (ielstart_all(ii-1,iproc).le.jj+1.and.
1279      &        ielend_all(ii-1,iproc).ge.jj+1) then
1280             flag=.true.
1281             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1282      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1283               isent=isent+1
1284               itask(isent)=iproc
1285               call add_task(iproc,ntask_cont_to,itask_cont_to)
1286             endif
1287           endif
1288         endif
1289       enddo
1290       return
1291       end
1292 c---------------------------------------------------------------------------
1293       subroutine add_int_from(ii,jj,ntask_cont_from,itask_cont_from)
1294       implicit none
1295       include "DIMENSIONS"
1296       include "COMMON.INTERACT"
1297       include "COMMON.SETUP"
1298       include "COMMON.IOUNITS"
1299       integer ii,jj,itask(2),ntask_cont_from,
1300      & itask_cont_from(0:max_fg_procs-1)
1301       logical flag
1302       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1303      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1304       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1305      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1306      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1307      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1308      & ielend_all(maxres,0:max_fg_procs-1)
1309       integer iproc,k,l
1310       do iproc=fg_rank+1,nfgtasks-1
1311         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1312           l=k+2
1313           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1 
1314      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1) 
1315      &    then
1316 c            write (iout,*)"turn3 from iproc",iproc," ij",ii,jj," kl",k,l
1317             call add_task(iproc,ntask_cont_from,itask_cont_from)
1318           endif
1319         enddo 
1320         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1321           l=k+3
1322           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1 
1323      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1) 
1324      &    then
1325 c            write (iout,*)"turn4 from iproc",iproc," ij",ii,jj," kl",k,l
1326             call add_task(iproc,ntask_cont_from,itask_cont_from)
1327           endif
1328         enddo 
1329         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0) then
1330           if (ii+1.ge.iatel_s_all(iproc).and.ii+1.le.iatel_e_all(iproc))
1331      &    then
1332             if (jj+1.ge.ielstart_all(ii+1,iproc).and.
1333      &          jj+1.le.ielend_all(ii+1,iproc)) then
1334               call add_task(iproc,ntask_cont_from,itask_cont_from)
1335             endif            
1336             if (jj-1.ge.ielstart_all(ii+1,iproc).and.
1337      &          jj-1.le.ielend_all(ii+1,iproc)) then
1338               call add_task(iproc,ntask_cont_from,itask_cont_from)
1339             endif
1340           endif
1341           if (ii-1.ge.iatel_s_all(iproc).and.ii-1.le.iatel_e_all(iproc))
1342      &    then
1343             if (jj-1.ge.ielstart_all(ii-1,iproc).and.
1344      &          jj-1.le.ielend_all(ii-1,iproc)) then
1345               call add_task(iproc,ntask_cont_from,itask_cont_from)
1346             endif
1347             if (jj+1.ge.ielstart_all(ii-1,iproc).and.
1348      &          jj+1.le.ielend_all(ii-1,iproc)) then
1349                call add_task(iproc,ntask_cont_from,itask_cont_from)
1350             endif
1351           endif
1352         endif
1353       enddo
1354       return
1355       end
1356 c---------------------------------------------------------------------------
1357       subroutine add_task(iproc,ntask_cont,itask_cont)
1358       implicit none
1359       include "DIMENSIONS"
1360       integer iproc,ntask_cont,itask_cont(0:max_fg_procs-1)
1361       integer ii
1362       do ii=1,ntask_cont
1363         if (itask_cont(ii).eq.iproc) return
1364       enddo
1365       ntask_cont=ntask_cont+1
1366       itask_cont(ntask_cont)=iproc
1367       return
1368       end
1369 c---------------------------------------------------------------------------
1370       subroutine int_bounds(total_ints,lower_bound,upper_bound)
1371       implicit real*8 (a-h,o-z)
1372       include 'DIMENSIONS'
1373       include 'mpif.h'
1374       include 'COMMON.SETUP'
1375       integer total_ints,lower_bound,upper_bound
1376       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1377       nint=total_ints/nfgtasks
1378       do i=1,nfgtasks
1379         int4proc(i-1)=nint
1380       enddo
1381       nexcess=total_ints-nint*nfgtasks
1382       do i=1,nexcess
1383         int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1384       enddo
1385       lower_bound=0
1386       do i=0,fg_rank-1
1387         lower_bound=lower_bound+int4proc(i)
1388       enddo 
1389       upper_bound=lower_bound+int4proc(fg_rank)
1390       lower_bound=lower_bound+1
1391       return
1392       end
1393 c---------------------------------------------------------------------------
1394       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1395       implicit real*8 (a-h,o-z)
1396       include 'DIMENSIONS'
1397       include 'mpif.h'
1398       include 'COMMON.SETUP'
1399       integer total_ints,lower_bound,upper_bound
1400       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1401       nint=total_ints/nfgtasks1
1402       do i=1,nfgtasks1
1403         int4proc(i-1)=nint
1404       enddo
1405       nexcess=total_ints-nint*nfgtasks1
1406       do i=1,nexcess
1407         int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1408       enddo
1409       lower_bound=0
1410       do i=0,fg_rank1-1
1411         lower_bound=lower_bound+int4proc(i)
1412       enddo 
1413       upper_bound=lower_bound+int4proc(fg_rank1)
1414       lower_bound=lower_bound+1
1415       return
1416       end
1417 c---------------------------------------------------------------------------
1418       subroutine int_partition(int_index,lower_index,upper_index,atom,
1419      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
1420       implicit real*8 (a-h,o-z)
1421       include 'DIMENSIONS'
1422       include 'COMMON.IOUNITS'
1423       integer int_index,lower_index,upper_index,atom,at_start,at_end,
1424      & first_atom,last_atom,int_gr,jat_start,jat_end
1425       logical lprn
1426       lprn=.false.
1427       if (lprn) write (iout,*) 'int_index=',int_index
1428       int_index_old=int_index
1429       int_index=int_index+last_atom-first_atom+1
1430       if (lprn) 
1431      &   write (iout,*) 'int_index=',int_index,
1432      &               ' int_index_old',int_index_old,
1433      &               ' lower_index=',lower_index,
1434      &               ' upper_index=',upper_index,
1435      &               ' atom=',atom,' first_atom=',first_atom,
1436      &               ' last_atom=',last_atom
1437       if (int_index.ge.lower_index) then
1438         int_gr=int_gr+1
1439         if (at_start.eq.0) then
1440           at_start=atom
1441           jat_start=first_atom-1+lower_index-int_index_old
1442         else
1443           jat_start=first_atom
1444         endif
1445         if (lprn) write (iout,*) 'jat_start',jat_start
1446         if (int_index.ge.upper_index) then
1447           at_end=atom
1448           jat_end=first_atom-1+upper_index-int_index_old
1449           return1
1450         else
1451           jat_end=last_atom
1452         endif
1453         if (lprn) write (iout,*) 'jat_end',jat_end
1454       endif
1455       return
1456       end
1457 #endif
1458 c------------------------------------------------------------------------------
1459       subroutine hpb_partition
1460       implicit real*8 (a-h,o-z)
1461       include 'DIMENSIONS'
1462 #ifdef MPI
1463       include 'mpif.h'
1464 #endif
1465       include 'COMMON.SBRIDGE'
1466       include 'COMMON.IOUNITS'
1467       include 'COMMON.SETUP'
1468 #ifdef MPI
1469       call int_bounds(nhpb,link_start,link_end)
1470       write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1471      &  ' absolute rank',MyRank,
1472      &  ' nhpb',nhpb,' link_start=',link_start,
1473      &  ' link_end',link_end
1474 #else
1475       link_start=1
1476       link_end=nhpb
1477 #endif
1478       return
1479       end