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