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