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