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