source/maxlik/src-Fmatch_safe/readrtns_forcefit.F

   1       subroutine atom_partition(iprot)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       include "COMMON.IOUNITS"
   6       include "COMMON.FMATCH"
   7       include "COMMON.INTERACT"
   8       include "COMMON.NAMES"
   9       include "COMMON.CHAIN"
  10       include "COMMON.PROTFILES"
  11       character*4 recid,atnam,resnam
  12       integer iprot,ires_prev,ires,nnres,nsc,iat,i,ii,j
  13       integer ilen
  14       external ilen
  15       logical carflag
  16       write (iout,*) "In atom_partition"
  17       call restore_molinfo(iprot)
  18       write (iout,*) "nres",nres
  19       open(ipdbin,file=fpdbfile(iprot),status="old",err=10)
  20       goto 11
  21    10 write (iout,*) "Error opening all-atom-type file ",
  22      &  fpdbfile(iprot)(:ilen(fpdbfile(iprot)))
  23       stop
  24    11 continue
  25       write (iout,*) "atom_partitio: after opening file"
  26 #ifdef PISIORNIA
  27       do i=1,nres-1
  28         print *,"i",i
  29         read(ipdbin,'(i3,5x,i3,1x,$)',err=20,end=20) iii,natp(i,iprot)
  30         do j=1,natp(i,iprot)
  31           read(ipdbin,'(a4,i4,1x,$)',err=20,end=20) atp(j,i,iprot),
  32      &         iatp(j,i,iprot)
  33         enddo
  34         read(ipdbin,'(3x,i3,1x,$)',err=20,end=20) natsc(i,iprot)
  35         do j=1,natsc(i,iprot)
  36           read(ipdbin,'(a4,i4,1x,$)',err=20,end=20) atsc(j,i,iprot),
  37      &         iatsc(j,i,iprot)
  38         enddo
  39         read(ipdbin,*)
  40       enddo
  41 #else
  42       natp(:,iprot)=0
  43       natsc(:,iprot)=1
  44       natsc(1,iprot)=0
  45       natsc(nres,iprot)=0
  46       ires_prev=0
  47       nat(iprot)=0
  48       nsc=0
  49       do
  50         read (ipdbin,'(a4,i7,1x,a4,1x,a3,i6)',end=20,err=20)
  51      &      recid,iat,atnam,resnam,ires
  52         write (iout,*) recid,iat,atnam,resnam,ires
  53         if (recid.ne."ATOM") cycle
  54         if (iat.gt.nat(iprot)) nat(iprot)=iat
  55         if (ires.ne.ires_prev) then
  56           carflag=.false.
  57           ires_prev=ires
  58         endif
  59         if (atnam.eq." CA ") then
  60           if (restyp(itype(ires+nnt-1)).ne.resnam) then
  61             write (iout,*)
  62      &     "Error: different residue name in force-match template",
  63      &     " protein",iprot," residue",ires,resnam,restyp(itype(ires+1))
  64             stop
  65           endif
  66           carflag=.true.
  67           iatsc(1,ires+nnt-1,iprot)=iat
  68           atsc(1,ires+nnt-1,iprot)=atnam(2:)
  69         endif
  70 c        print *,"carflag ",carflag
  71         if (.not.carflag) then
  72           if (resnam.eq."PRO ") then
  73             write (iout,*) resnam,ires
  74 c            if (atnam.eq." N  ".or.atnam.eq." CD ".or.atnam.eq." HD2"
  75 c            if (atnam.eq." N  ".or.atnam.eq." HD2".or.atnam.eq." HD3")
  76             if (atnam.eq." N  ")
  77      &      then
  78               natp(ires-1+nnt-1,iprot)=natp(ires-1+nnt-1,iprot)+1
  79 c              print *,"back",natp(ires-1)
  80               ii=natp(ires-1+nnt-1,iprot)
  81               atp(ii,ires-1+nnt-1,iprot)=atnam(2:)
  82               iatp(ii,ires-1+nnt-1,iprot)=iat
  83             else
  84               natsc(ires+nnt-1,iprot)=natsc(ires+nnt-1,iprot)+1
  85 c              print *,"sc",natsc(ires)
  86               ii=natsc(ires+nnt-1,iprot)
  87               iatsc(ii,ires+nnt-1,iprot)=iat
  88               atsc(ii,ires+nnt-1,iprot)=atnam(2:)
  89             endif
  90           else
  91             natp(ires-1+nnt-1,iprot)=natp(ires-1+nnt-1,iprot)+1
  92             ii=natp(ires-1+nnt-1,iprot)
  93             atp(ii,ires-1+nnt-1,iprot)=atnam(2:)
  94             iatp(ii,ires-1+nnt-1,iprot)=iat
  95           endif
  96         elseif(atnam.eq." C  ".or.atnam.eq." O  " .or.atnam.eq." OXT")
  97      &  then
  98           natp(ires+nnt-1,iprot)=natp(ires+nnt-1,iprot)+1
  99           ii=natp(ires+nnt-1,iprot)
 100           if (atnam(1:1).eq." ") then
 101             atp(ii,ires+nnt-1,iprot)=atnam(2:)
 102           else
 103             atp(ii,ires+nnt-1,iprot)=atnam
 104           endif
 105           iatp(ii,ires+nnt-1,iprot)=iat
 106         else if (atnam.ne. " CA ") then
 107           natsc(ires+nnt-1,iprot)=natsc(ires+nnt-1,iprot)+1
 108           ii=natsc(ires+nnt-1,iprot)
 109           iatsc(ii,ires+nnt-1,iprot)=iat
 110           if (atnam(1:1).eq." ") then
 111             atsc(ii,ires+nnt-1,iprot)=atnam(2:)
 112           else
 113             atsc(ii,ires+nnt-1,iprot)=atnam
 114           endif
 115         endif
 116         nnres=ires
 117         if (atnam.eq." CA ".and.resnam.ne."GLY ") nsc=nsc+1
 118       enddo
 119    20 continue
 120       close(ipdbin)
 121       write (iout,*) "End reading nnres",nnres," nres",nres,
 122      &   " nsc",nsc," nat",nat(iprot)
 123 #endif
 124       close(ipdbin)
 125 c---------------------------------------------------------------------
 126       do i=1,nres
 127         write(iout,'(i3,1x,a3,1x,5h back,i3,1x,$)')i,restyp(itype((i))),
 128      &   natp(i,iprot)
 129         do j=1,natp(i,iprot)
 130           write(iout,'(a4,i4,1x,$)') atp(j,i,iprot),iatp(j,i,iprot)
 131         enddo
 132         write(iout,'(i3,1x,3h sc,i3,1x,$)') i,natsc(i,iprot)
 133         do j=1,natsc(i,iprot)
 134           write(iout,'(a4,i4,1x,$)') atsc(j,i,iprot),iatsc(j,i,iprot)
 135         enddo
 136         write (iout,*)
 137       enddo
 138       nsite(iprot)=nres
 139       do i=nnt,nct
 140         if (itype(i).ne.10) nsite(iprot)=nsite(iprot)+1
 141       enddo
 142       write (iout,*) "iprot",iprot," nsite",nsite(iprot),nres+nsc
 143       return
 144       end
 145 c---------------------------------------------------------------------
 146       subroutine read_allat_coord_forces(iprot)
 147       implicit none
 148       include "DIMENSIONS"
 149       include "DIMENSIONS.ZSCOPT"
 150 #ifdef MPI
 151       include "mpif.h"
 152       include "COMMON.MPI"
 153       integer ierror,errcode
 154 #endif
 155       include "COMMON.CONTROL"
 156       include "COMMON.IOUNITS"
 157       include "COMMON.FMATCH"
 158       include "COMMON.CHAIN"
 159       include "COMMON.INTERACT"
 160       include "COMMON.PROTFILES"
 161       include "COMMON.SBRIDGE"
 162       include "COMMON.NAMES"
 163       include "COMMON.WEIGHTS"
 164       include "COMMON.CLASSES"
 165       include "COMMON.PROTNAME"
 166       integer i,ii,j,k,l,kkk,kkkk,jj_old,jj,jjj,icount,nnsc,iis,jjs
 167       integer iprot
 168       integer ilen
 169       external ilen
 170       double precision casc(3),dss
 171       real rjunk
 172
 173       jj_old=1
 174       call restore_molinfo(iprot)
 175       open (ientout,file=bprotfiles_MD(iprot),status="unknown",
 176      &    form="unformatted",access="direct",recl=lenrec_MD(iprot))
 177 c Change AL 12/30/2017
 178 #ifdef MPI
 179       if (me.eq.Master) then
 180 ! Loop over snapshots
 181 #endif!
 182       open(ipdbin,file=fcoordfile(iprot),status="old",err=10)
 183       goto 11
 184    10 write (iout,*) "Error opening MD coordinate file ",
 185      &  fcoordfile(iprot)(:ilen(fcoordfile(iprot)))
 186       call flush(iout)
 187 #ifdef MPI
 188       call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR)
 189 #else
 190       stop
 191 #endif
 192    11 continue
 193       open(ientin,file=fforcefile(iprot),status="old",err=20)
 194       goto 21
 195    20 write (iout,*) "Error opening MD force file ",
 196      &  fforcefile(iprot)(:ilen(fforcefile(iprot)))
 197       call flush(iout)
 198 #ifdef MPI
 199       call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR)
 200 #else
 201       stop
 202 #endif
 203    21 continue
 204       write (iout,*) "nat",nat(iprot)
 205       call flush(iout)
 206       kkkk=0
 207       jjj=0
 208       jj=0
 209       icount=0
 210       CGforce_MD_ave(:,:,:,iprot)=0.0d0
 211       write (iout,*) "istart_MD",istart_MD(iprot),
 212      &   " iend_MD",iend_MD(iprot)
 213       read (ipdbin,*)
 214       read (ientin,*)
 215       do kkk=1,iend_MD(iprot)
 216 c      do kkk=1,1
 217       read (ipdbin,'(10f8.3)',end=30,err=30)
 218      &      ((atcoord(j,i),j=1,3),i=1,nat(iprot))
 219       if (mdbox(iprot)) read(ipdbin,*,end=30,err=30) rjunk,rjunk,rjunk
 220       read(ientin,'(10f8.3)',end=30,err=30)((atforce(j,i),j=1,3),i=1,
 221      &      nat(iprot))
 222       if (kkk.lt.istart_MD(iprot)) cycle
 223       call compute_CG_forces_ave(kkkk+1,iprot)
 224       if (mod(kkk,ifreq_MD(iprot)).gt.0) cycle
 225       kkkk = kkkk + 1
 226       do i=1,2*nres
 227         do j=1,3
 228           CGforce_MD_ave(j,i,kkkk,iprot)=CGforce_MD_ave(j,i,kkkk,iprot)/
 229      &      ifreq_MD(iprot)
 230         enddo
 231       enddo
 232 #ifdef DEBUG
 233       write(iout,*) "All-atom Coordinates and forces",kkk,kkkk
 234       do i=1,nres-1
 235         do j=1,natp(i,iprot)
 236           write(iout,'(2i3,1x,a4,i5,3f8.3,5x,3f8.3)')
 237      &     i,j,atp(j,i,iprot),iatp(j,i,iprot),
 238      &     (atcoord(k,iatp(j,i,iprot)),k=1,3),
 239      &     (atforce(k,iatp(j,i,iprot)),k=1,3)
 240         enddo
 241         do j=1,natsc(i,iprot)
 242          write(iout,'(2i3,1x,a4,i5,3f8.3,5x,3f8.3)')
 243      &     i,j,atsc(j,i,iprot),iatsc(j,i,iprot),
 244      &     (atcoord(k,iatsc(j,i,iprot)),k=1,3),
 245      &     (atforce(k,iatsc(j,i,iprot)),k=1,3)
 246         enddo
 247       enddo
 248 #endif
 249       do k=1,3
 250         c(k,1)=atcoord(k,iatp(1,1,iprot))
 251         c(k,nres)=atcoord(k,iatp(3,nres-1,iprot))
 252       enddo
 253       do i=nnt,nct
 254         do j=1,3
 255           c(j,i)=atcoord(j,iatsc(1,i,iprot))
 256         enddo
 257       enddo
 258       iis=0
 259       do i=1,nres
 260         casc=0.0d0
 261 c        write (iout,*) "i",i," natsc",natsc(i)
 262         nnsc=0
 263         do j=1,natsc(i,iprot)
 264 c          write (iout,*) "j",j," iatsc",iatsc(j,i)
 265           if (atsc(j,i,iprot)(1:1).eq."H") cycle
 266           nnsc=nnsc+1
 267           do k=1,3
 268             casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot))
 269           enddo
 270         enddo
 271         do k=1,3
 272           c(k,i+nres)=casc(k)/nnsc
 273         enddo
 274 c        write (iout,*) (c(k,i+nres),k=1,3)
 275       enddo
 276 c Detect disulfide bonds
 277       nss=0
 278       do i=nnt,nct
 279         do j=i+1,nct
 280           if (itype(i).eq.1 .and. itype(j).eq.1) then
 281             do k=1,natsc(i,iprot)
 282               do l=1,natsc(j,iprot)
 283                 if (atsc(k,i,iprot).eq." S  ".and.
 284      &              atsc(l,j,iprot).eq." S  ") then
 285                   iis=iatsc(k,i,iprot)
 286                   jjs=iatsc(l,j,iprot)
 287                   dSS=dsqrt((atcoord(1,iis)-atcoord(1,jjs))**2+
 288      &                      (atcoord(2,iis)-atcoord(2,jjs))**2+
 289      &                      (atcoord(3,iis)-atcoord(3,jjs))**2)
 290                 endif
 291                 if (dSS.lt.2.5d0) then
 292                   nss=nss+1
 293                   ihpb(nss)=iis
 294                   jhpb(nss)=jjs
 295                 endif
 296               enddo
 297             enddo
 298           endif
 299         enddo
 300       enddo
 301 #ifdef DEBUG
 302       write (iout,*) "CG Coordinates"
 303       do i=1,nres
 304         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
 305           write (iout,'(a4,i5,3f10.3)')restyp(itype(i)),i,(c(j,i),j=1,3)
 306         else
 307           write (iout,'(a4,i5,3f10.3,5x,3f10.3)')
 308      &       restyp(itype(i)),i,(c(j,i),j=1,3),(c(j,i+nres),j=1,3)
 309         endif
 310       enddo
 311       write (iout,'(8f10.5)')
 312      &       ((c(l,k),l=1,3),k=1,nres),
 313      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 314       call int_from_cart1(.true.)
 315 #endif
 316
 317       call add_new_MDconf(jjj,jj,jj_old,icount,iprot)
 318 #ifdef DEBUG
 319       write (iout,*) "Average and momentary CG forces",
 320      &  " calculated from all-atom forces"
 321       write (iout,"(4x,a,t55,a)") "Backbone","Sidechain"
 322       ii=nres
 323       do i=1,nres
 324         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
 325           write (iout,'(a4,i5,3e15.5)') restyp(itype(i)),i,
 326      &      (CGforce_MD_ave(j,i,kkkk,iprot),j=1,3)
 327           write (iout,'(9x,3e15.5)')
 328      &      (CGforce_MD(j,i,kkkk,iprot),j=1,3)
 329         else
 330           ii=ii+1
 331           write (iout,'(a4,i5,3e15.5,5x,3e15.5)') restyp(itype(i)),i,
 332      &      (CGforce_MD_ave(j,i,kkkk,iprot),j=1,3),
 333      &      (CGforce_MD_ave(j,ii,kkkk,iprot),j=1,3)
 334           write (iout,'(9x,3e15.5,5x,3e15.5)')
 335      &      (CGforce_MD(j,i,kkkk,iprot),j=1,3),
 336      &      (CGforce_MD(j,ii,kkkk,iprot),j=1,3)
 337         endif
 338       enddo
 339 #endif
 340
 341       enddo
 342    30 continue
 343       if (icount.gt.0) then
 344       call write_and_send_MDconf(icount,jj_old,jj,iprot)
 345       nconf_forc(iprot)=kkkk
 346       ntot_work_MD(iprot)=kkkk
 347       endif
 348       call write_and_send_MDconf(0,jj_old,jj,iprot)
 349 #ifdef MPI
 350       else
 351 #endif
 352       jjj=0
 353       jj=0
 354       icount=0
 355       do i=1,maxstr
 356         list_conf_MD(i,iprot)=i
 357       enddo
 358       call receive_and_write_MDconf(icount,jj_old,jj,iprot)
 359 #ifdef MPI
 360       endif
 361       close(icbase)
 362       ntot_MD(iprot)=jj
 363         write (iout,*) "Protein",
 364      &    protname(iprot)(:ilen(protname(iprot))),", ",ntot_MD(iprot),
 365      &   " conformatons read ",ntot_MD(iprot),
 366      &   " conformations written to ",
 367      &  bprotfiles_MD(iprot)(:ilen(bprotfiles_MD(iprot)))
 368         ntot_MD(0) = ntot_MD(0)+ntot_MD(iprot)
 369 #endif
 370       return
 371       end
 372 c-------------------------------------------------------------------
 373       subroutine compute_CG_forces(kkkk,iprot)
 374       implicit none
 375       include "DIMENSIONS"
 376       include "DIMENSIONS.ZSCOPT"
 377       include "COMMON.IOUNITS"
 378       include "COMMON.FMATCH"
 379       include "COMMON.CHAIN"
 380       include "COMMON.INTERACT"
 381       include "COMMON.NAMES"
 382       integer i,ii,iis,j,k,kkkk,iprot,nsc,nscc,nnsc,n,nn
 383       double precision yy(3),caca(3),casc(3),xx,cacanorm,cascnorm
 384       double precision scalar
 385       ii=0
 386       CGforce_MD(:,:,kkkk,iprot) = 0.0d0
 387       do j=1,natp(1,iprot)
 388 c        write (iout,'(3i5,3f10.5)') kkkk,j,iatp(j,1,iprot),
 389 c     &     (atforce(k,iatp(j,1,iprot)),k=1,3)
 390         do k=1,3
 391           CGforce_MD(k,nnt,kkkk,iprot)=CGforce_MD(k,nnt,kkkk,iprot)
 392      &      +atforce(k,iatp(j,1,iprot))
 393         enddo
 394       enddo
 395       do i=nnt,nct-1
 396         do j=1,3
 397           caca(j)=atcoord(j,iatsc(1,i+1,iprot))
 398      &        -atcoord(j,iatsc(1,i,iprot))
 399         enddo
 400         cacanorm=dsqrt(scalar(caca,caca))
 401 c        write(iout,*)"i",i," cacanorm",cacanorm
 402         do k=1,3
 403           caca(k)=caca(k)/cacanorm
 404         enddo
 405         do j=1,natp(i,iprot)
 406           yy=0.0d0
 407           do k=1,3
 408             yy(k)=atcoord(k,iatp(j,i,iprot))-atcoord(k,iatsc(1,i,iprot))
 409           enddo
 410           xx=scalar(caca,yy)/cacanorm
 411           do k=1,3
 412             CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot)
 413      &        +(1-xx)*atforce(k,iatp(j,i,iprot))
 414             CGforce_MD(k,i+1,kkkk,iprot)=CGforce_MD(k,i+1,kkkk,iprot)
 415      &        +xx*atforce(k,iatp(j,i,iprot))
 416           enddo
 417         enddo
 418       enddo
 419 c#ifdef SC
 420       nscc=0
 421       iis=nres
 422       do i=nnt,nct
 423         casc=0.0d0
 424 c        print *,"i",i," natsc",natsc(i)
 425         nnsc=0
 426         do j=1,natsc(i,iprot)
 427 c          print *,"j",j," iatsc",iatsc(j,i)
 428           if (atsc(j,i,iprot)(1:1).eq."H") cycle
 429           nnsc=nnsc+1
 430           do k=1,3
 431             casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot))
 432           enddo
 433         enddo
 434         do k=1,3
 435           casc(k)=casc(k)/nnsc-atcoord(k,iatsc(1,i,iprot))
 436         enddo
 437         cascnorm=dsqrt(scalar(casc,casc))
 438 c        write(iout,*)"i",i,restyp(itype(i))," natsc",natsc(i,iprot),
 439 c     &    " nnsc",nnsc," cascnorm",cascnorm
 440         if (nnsc.gt.1) then
 441           nscc=nscc+1
 442           iis=iis+1
 443           ires_sc(iis,iprot)=i
 444         endif
 445 c        write (iout,*) "nnsc",nnsc," iis",iis
 446         if (cascnorm.gt.1.0d-2) then
 447           do k=1,3
 448             casc(k)=casc(k)/cascnorm
 449           enddo
 450         else
 451           casc=0.0d0
 452         endif
 453         do k=1,3
 454           CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot)
 455      &      +atforce(k,iatsc(1,i,iprot))
 456         enddo
 457         do j=2,natsc(i,iprot)
 458
 459           if (atsc(j,i,iprot)(1:2).eq."HA") then
 460
 461           do k=1,3
 462             CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot)
 463      &       +atforce(k,iatsc(j,i,iprot))
 464           enddo
 465 c          write (4,*) "HA",i,(CGforce_MD(k,i,kkkk),k=1,3)
 466
 467           else
 468
 469           do k=1,3
 470            yy(k)=atcoord(k,iatsc(j,i,iprot))-atcoord(k,iatsc(1,i,iprot))
 471           enddo
 472           xx=scalar(casc,yy)/cascnorm
 473           do k=1,3
 474             CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot)
 475      &       +(1.0d0-xx)*atforce(k,iatsc(j,i,iprot))
 476             CGforce_MD(k,iis,kkkk,iprot)=CGforce_MD(k,iis,kkkk,iprot)
 477      &       +xx*atforce(k,iatsc(j,i,iprot))
 478           enddo
 479 c          write (4,*)"REG CA",i,(CGforce_MD(k,i,kkkk),k=1,3)
 480 c          write (4,*)"REG SC",iis,(CGforce_MD(k,iis,kkkk),k=1,3)
 481
 482           endif
 483
 484         enddo
 485       enddo
 486       do j=1,natp(nct,iprot)
 487         do k=1,3
 488           CGforce_MD(k,nct,kkkk,iprot)=CGforce_MD(k,nct,kkkk,iprot)
 489      &     +atforce(k,iatp(j,nct,iprot))
 490         enddo
 491       enddo
 492 c#endif
 493       nn = ii
 494 c      print *,"nres",nres," nsc",nsc
 495       n = nres+nsc
 496 #ifdef DEBUG
 497       write(iout,*)"CG forces calculated from all-atom forces kkkk",kkkk
 498       write (iout,"(4x,a,t55,a)") "Backbone","Sidechain"
 499       ii=nres
 500       do i=1,nres
 501         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
 502           write (iout,'(a4,i5,3e15.5)') restyp(itype(i)),i,
 503      &      (CGforce_MD(j,i,kkkk,iprot),j=1,3)
 504         else
 505           ii=ii+1
 506           write (iout,'(a4,i5,3e15.5,5x,3e15.5)') restyp(itype(i)),i,
 507      &      (CGforce_MD(j,i,kkkk,iprot),j=1,3),
 508      &      (CGforce_MD(j,ii,kkkk,iprot),j=1,3)
 509         endif
 510       enddo
 511 #endif
 512       return
 513       end
 514 c-------------------------------------------------------------------
 515       subroutine compute_CG_forces_ave(kkkk,iprot)
 516       implicit none
 517       include "DIMENSIONS"
 518       include "DIMENSIONS.ZSCOPT"
 519       include "COMMON.IOUNITS"
 520       include "COMMON.FMATCH"
 521       include "COMMON.CHAIN"
 522       include "COMMON.INTERACT"
 523       include "COMMON.NAMES"
 524       integer i,ii,iis,j,k,kkkk,iprot,nsc,nscc,nnsc,n,nn
 525       double precision yy(3),caca(3),casc(3),xx,cacanorm,cascnorm
 526       double precision scalar
 527       ii=0
 528
 529       do j=1,natp(1,iprot)
 530 c        write (iout,'(3i5,3f10.5)') kkkk,j,iatp(j,1,iprot),
 531 c     &     (atforce(k,iatp(j,1,iprot)),k=1,3)
 532         do k=1,3
 533           CGforce_MD_ave(k,nnt,kkkk,iprot)=
 534      &      CGforce_MD_ave(k,nnt,kkkk,iprot)
 535      &      +atforce(k,iatp(j,1,iprot))
 536         enddo
 537       enddo
 538       do i=nnt,nct-1
 539         do j=1,3
 540           caca(j)=atcoord(j,iatsc(1,i+1,iprot))
 541      &        -atcoord(j,iatsc(1,i,iprot))
 542         enddo
 543         cacanorm=dsqrt(scalar(caca,caca))
 544 c        write(iout,*)"i",i," cacanorm",cacanorm
 545         do k=1,3
 546           caca(k)=caca(k)/cacanorm
 547         enddo
 548         do j=1,natp(i,iprot)
 549           yy=0.0d0
 550           do k=1,3
 551             yy(k)=atcoord(k,iatp(j,i,iprot))-atcoord(k,iatsc(1,i,iprot))
 552           enddo
 553           xx=scalar(caca,yy)/cacanorm
 554           do k=1,3
 555            CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD_ave(k,i,kkkk,iprot)
 556      &        +(1-xx)*atforce(k,iatp(j,i,iprot))
 557            CGforce_MD_ave(k,i+1,kkkk,iprot)=
 558      &        CGforce_MD_ave(k,i+1,kkkk,iprot)
 559      &        +xx*atforce(k,iatp(j,i,iprot))
 560           enddo
 561         enddo
 562       enddo
 563 c#ifdef SC
 564       nscc=0
 565       iis=nres
 566       do i=nnt,nct
 567         casc=0.0d0
 568 c        print *,"i",i," natsc",natsc(i)
 569         nnsc=0
 570         do j=1,natsc(i,iprot)
 571 c          print *,"j",j," iatsc",iatsc(j,i)
 572           if (atsc(j,i,iprot)(1:1).eq."H") cycle
 573           nnsc=nnsc+1
 574           do k=1,3
 575             casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot))
 576           enddo
 577         enddo
 578         do k=1,3
 579           casc(k)=casc(k)/nnsc-atcoord(k,iatsc(1,i,iprot))
 580         enddo
 581         cascnorm=dsqrt(scalar(casc,casc))
 582 c        write(iout,*)"i",i,restyp(itype(i))," natsc",natsc(i,iprot),
 583 c     &    " nnsc",nnsc," cascnorm",cascnorm
 584         if (nnsc.gt.1) then
 585           nscc=nscc+1
 586           iis=iis+1
 587           ires_sc(iis,iprot)=i
 588         endif
 589 c        write (iout,*) "nnsc",nnsc," iis",iis
 590         if (cascnorm.gt.1.0d-2) then
 591           do k=1,3
 592             casc(k)=casc(k)/cascnorm
 593           enddo
 594         else
 595           casc=0.0d0
 596         endif
 597         do k=1,3
 598           CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD_ave(k,i,kkkk,iprot)
 599      &      +atforce(k,iatsc(1,i,iprot))
 600         enddo
 601         do j=2,natsc(i,iprot)
 602
 603           if (atsc(j,i,iprot)(1:2).eq."HA") then
 604
 605           do k=1,3
 606            CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD_ave(k,i,kkkk,iprot)
 607      &       +atforce(k,iatsc(j,i,iprot))
 608           enddo
 609 c          write (4,*) "HA",i,(CGforce_MD(k,i,kkkk),k=1,3)
 610
 611           else
 612
 613           do k=1,3
 614            yy(k)=atcoord(k,iatsc(j,i,iprot))-atcoord(k,iatsc(1,i,iprot))
 615           enddo
 616           xx=scalar(casc,yy)/cascnorm
 617           do k=1,3
 618            CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD_ave(k,i,kkkk,iprot)
 619      &       +(1.0d0-xx)*atforce(k,iatsc(j,i,iprot))
 620            CGforce_MD_ave(k,iis,kkkk,iprot)=
 621      &       CGforce_MD_ave(k,iis,kkkk,iprot)
 622      &       +xx*atforce(k,iatsc(j,i,iprot))
 623           enddo
 624 c          write (4,*)"REG CA",i,(CGforce_MD(k,i,kkkk),k=1,3)
 625 c          write (4,*)"REG SC",iis,(CGforce_MD(k,iis,kkkk),k=1,3)
 626
 627           endif
 628
 629         enddo
 630       enddo
 631       do j=1,natp(nct,iprot)
 632         do k=1,3
 633           CGforce_MD_ave(k,nct,kkkk,iprot)=
 634      &     CGforce_MD_ave(k,nct,kkkk,iprot)
 635      &     +atforce(k,iatp(j,nct,iprot))
 636         enddo
 637       enddo
 638 c#endif
 639       nn = ii
 640 c      print *,"nres",nres," nsc",nsc
 641       n = nres+nsc
 642       return
 643       end
 644 c----------------------------------------------------------------------
 645       subroutine add_new_MDconf(jjj,jj,jj_old,icount,iprot)
 646       implicit none
 647       include "DIMENSIONS"
 648       include "DIMENSIONS.ZSCOPT"
 649       include "COMMON.CHAIN"
 650       include "COMMON.INTERACT"
 651       include "COMMON.LOCAL"
 652       include "COMMON.CLASSES"
 653       include "COMMON.ENERGIES"
 654       include "COMMON.IOUNITS"
 655       include "COMMON.PROTFILES"
 656       include "COMMON.PROTNAME"
 657       include "COMMON.VMCPAR"
 658       include "COMMON.WEIGHTS"
 659       include "COMMON.NAMES"
 660       include "COMMON.ALLPROT"
 661       include "COMMON.WEIGHTDER"
 662       include "COMMON.VAR"
 663       include "COMMON.SBRIDGE"
 664       include "COMMON.GEO"
 665       include "COMMON.COMPAR"
 666       integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,ib,ibatch,
 667      &  nn,nn1,inan,Next,itj
 668       jjj=jjj+1
 669 c      write(iout,*) "add_new_MDconf jjj jj jj_old:",jjj,jj,jj_old
 670 c      call flush(iout)
 671       call int_from_cart1(.false.)
 672 c      call flush(iout)
 673 c      write (iout,*) "After int_from_cart1"
 674       call flush(iout)
 675       do j=nnt+1,nct
 676         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.5d0) then
 677           write (iout,*) "nnt",nnt," nct",nct
 678           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j)
 679           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
 680           write (iout,*) "The Cartesian geometry is:"
 681           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 682           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 683           write (iout,*) "The internal geometry is:"
 684           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 685           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 686           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 687           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 688           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 689           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 690           write (iout,*)
 691      &      "This conformation WILL NOT be added to the database."
 692           return
 693         endif
 694       enddo
 695       do j=nnt,nct
 696         itj=itype(j)
 697         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
 698           write (iout,*) "nnt",nnt," nct",nct
 699           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
 700           write (iout,*) "The Cartesian geometry is:"
 701           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 702           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 703           write (iout,*) "The internal geometry is:"
 704           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 705           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 706           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 707           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 708           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 709           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 710           write (iout,*)
 711      &      "This conformation WILL NOT be added to the database."
 712           return
 713         endif
 714       enddo
 715       do j=3,nres
 716         if (theta(j).le.0.0d0) then
 717           write (iout,*)
 718      &      "Zero theta angle(s) in conformation",ii
 719           write (iout,*) "The Cartesian geometry is:"
 720           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 721           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 722           write (iout,*) "The internal geometry is:"
 723           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 724           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 725           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 726           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 727           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 728           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 729           write (iout,*)
 730      &      "This conformation WILL NOT be added to the database."
 731           return
 732         endif
 733         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 734       enddo
 735       jj=jj+1
 736       icount=icount+1
 737       list_conf_MD(jj,iprot)=jj
 738 c      write (iout,*) "calling store_MDconf",jjj,jj,icount
 739 c      call flush(iout)
 740       call store_MDconf_from_file(jj,icount,iprot)
 741 c      write (iout,*) "finished store_MDconf"
 742 c      call flush(iout)
 743       if (icount.eq.maxstr_proc) then
 744 #ifdef DEBUG
 745         write (iout,* ) "jj_old",jj_old," jj",jj
 746         write (iout,*) "Calling write_and_send_cconf"
 747         call flush(iout)
 748 #endif
 749         call write_and_send_MDconf(icount,jj_old,jj,iprot)
 750         jj_old=jj+1
 751 #ifdef DEBUG
 752         write (iout,*) "Exited write_and_send_cconf"
 753         call flush(iout)
 754 #endif
 755         icount=0
 756       endif
 757       return
 758       end
 759 c------------------------------------------------------------------------------
 760       subroutine store_MDconf_from_file(jj,icount,iprot)
 761       implicit none
 762       include "DIMENSIONS"
 763       include "DIMENSIONS.ZSCOPT"
 764       include "COMMON.CHAIN"
 765       include "COMMON.SBRIDGE"
 766       include "COMMON.INTERACT"
 767       include "COMMON.IOUNITS"
 768       include "COMMON.CLASSES"
 769       include "COMMON.ALLPROT"
 770       include "COMMON.VAR"
 771       include "COMMON.FMATCH"
 772       integer i,j,jj,icount,ibatch,iprot
 773 c Store the conformation that has been read in
 774       do i=1,2*nres
 775         do j=1,3
 776           cg_MD(j,i,icount,iprot)=c(j,i)
 777         enddo
 778       enddo
 779       nss_MD(icount,iprot)=nss
 780       do i=1,nss
 781         ihpb_MD(i,icount,iprot)=ihpb(i)
 782         jhpb_MD(i,icount,iprot)=jhpb(i)
 783       enddo
 784       call compute_CG_forces(icount,iprot)
 785       return
 786       end
 787 c------------------------------------------------------------------------------
 788       subroutine write_and_send_MDconf(icount,jj_old,jj,iprot)
 789       implicit none
 790       include "DIMENSIONS"
 791       include "DIMENSIONS.ZSCOPT"
 792 #ifdef MPI
 793       include "mpif.h"
 794       integer IERROR
 795       include "COMMON.MPI"
 796 #endif
 797       include "COMMON.WEIGHTS"
 798       include "COMMON.CHAIN"
 799       include "COMMON.SBRIDGE"
 800       include "COMMON.INTERACT"
 801       include "COMMON.IOUNITS"
 802       include "COMMON.CLASSES"
 803       include "COMMON.VAR"
 804       include "COMMON.ALLPROT"
 805       include "COMMON.ENERGIES"
 806       include "COMMON.WEIGHTDER"
 807       include "COMMON.OPTIM"
 808       include "COMMON.COMPAR"
 809       include "COMMON.FMATCH"
 810       integer i,k,l,icount,jj_old,jj,Next,iprot
 811 c Write the structures to a scratch file
 812 #ifdef MPI
 813 #ifdef DEBUG
 814       write (iout,*) "Processor",me," entered WRITE_AND_SEND_MDCONF"
 815       write (iout,*) "iprot",iprot
 816       call flush(iout)
 817 #endif
 818       call MPI_Bcast(icount,1,MPI_INTEGER,Master,WHAM_COMM,IERROR)
 819 #ifdef DEBUG
 820       write (iout,*) "Processor",me," sent icount=",icount
 821       write (iout,*) "Processor",me," jj_old",jj_old," jj",jj
 822       call flush(iout)
 823 #endif
 824       if (icount.gt.0) then
 825       call MPI_Bcast(jj_old,1,MPI_INTEGER,Master,WHAM_COMM,IERROR)
 826       call MPI_Bcast(jj,1,MPI_INTEGER,Master,WHAM_COMM,IERROR)
 827       write (iout,*) "Finished broadcast jj_old jj"
 828       call flush(iout)
 829       call MPI_Bcast(nss_MD(1,iprot),icount,MPI_INTEGER,Master,
 830      &    WHAM_COMM,IERROR)
 831       write (iout,*) "Finished broadcast nss"
 832       call flush(iout)
 833       call MPI_Bcast(ihpb_MD(1,1,iprot),maxres*icount,MPI_INTEGER,
 834      &    Master,WHAM_COMM,IERROR)
 835       call MPI_Bcast(jhpb_MD(1,1,iprot),maxres*icount,MPI_INTEGER,
 836      &    Master,WHAM_COMM,IERROR)
 837       write (iout,*) "Finished broadcast ihpb jhpb"
 838       call flush(iout)
 839       call MPI_Bcast(cg_MD(1,1,1,iprot),3*icount*maxres2,
 840      &    MPI_REAL,Master,WHAM_COMM,IERROR)
 841       write (iout,*) "Finished broadcast cg_MD"
 842       call flush(iout)
 843       call MPI_Bcast(cgforce_MD(1,1,1,iprot),3*icount*maxres2,
 844      &    MPI_REAL,Master,WHAM_COMM,IERROR)
 845       write (iout,*) "Finished broadcast cgforce_MD"
 846       call flush(iout)
 847       endif
 848       write (iout,*) "Master finished broadcast"
 849       call flush(iout)
 850 #endif
 851       call dawrite_MDcoords(iprot,jj_old,jj,ientout)
 852       nconf_forc(iprot)=jj
 853       write (iout,*) "iprot",iprot," jj",jj," nconf_forc",
 854      & nconf_forc(iprot)
 855 #ifdef DEBUG
 856       write (iout,*) "Protein",iprot
 857       write(iout,*)"Processor",me," received",icount,
 858      &  " MD conformations and forces"
 859       do i=1,icount
 860         write (iout,'(8f10.4)') ((cg_MD(l,k,i,iprot),l=1,3),k=1,nres)
 861         write(iout,'(8f10.4)')((cg_MD(l,k+nres,i,iprot),l=1,3),
 862      &    k=nnt,nct)
 863         write (iout,'(16i5)') nss_MD(i,iprot),(ihpb_MD(k,i,iprot),
 864      &    jhpb_MD(k,i,iprot),k=1,nss_MD(i,iprot))
 865         write(iout,'(8f10.4)')((cgforce_MD(l,k,i,iprot),l=1,3),k=1,
 866      &    nsite(iprot))
 867       enddo
 868 #endif
 869       return
 870       end
 871 c------------------------------------------------------------------------------
 872 #ifdef MPI
 873       subroutine receive_and_write_MDconf(icount,jj_old,jj,iprot)
 874       implicit none
 875       include "DIMENSIONS"
 876       include "DIMENSIONS.ZSCOPT"
 877       include "mpif.h"
 878       integer IERROR,STATUS(MPI_STATUS_SIZE)
 879       include "COMMON.MPI"
 880       include "COMMON.CHAIN"
 881       include "COMMON.SBRIDGE"
 882       include "COMMON.INTERACT"
 883       include "COMMON.IOUNITS"
 884       include "COMMON.CLASSES"
 885       include "COMMON.FMATCH"
 886       include "COMMON.ENERGIES"
 887       include "COMMON.VAR"
 888       include "COMMON.GEO"
 889       include "COMMON.WEIGHTS"
 890       include "COMMON.WEIGHTDER"
 891       include "COMMON.COMPAR"
 892       include "COMMON.OPTIM"
 893       integer i,j,k,l,icount,jj_old,jj,iprot,Previous,Next
 894       icount=1
 895       do while (icount.gt.0)
 896 #ifdef DEBUG
 897       write (iout,*) "Processor",me," entered RECEIVE_AND_SEND_MDCONF"
 898       write (iout,*) "iprot",iprot
 899       call flush(iout)
 900 #endif
 901       call MPI_Bcast(icount,1,MPI_INTEGER,Master,WHAM_COMM,IERROR)
 902 #ifdef DEBUG
 903       write (iout,*) "Processor",me," received icount=",icount
 904       call flush(iout)
 905 #endif
 906       if (icount.gt.0) then
 907       call MPI_Bcast(jj_old,1,MPI_INTEGER,Master,WHAM_COMM,IERROR)
 908       call MPI_Bcast(jj,1,MPI_INTEGER,Master,WHAM_COMM,IERROR)
 909       write (iout,*) "Finished broadcast jj_old jj"
 910       write (iout,*) "Processor",me," jj_old",jj_old," jj",jj
 911       call flush(iout)
 912       call MPI_Bcast(nss_MD(1,iprot),icount,MPI_INTEGER,Master,
 913      &    WHAM_COMM,IERROR)
 914       write (iout,*) "Finished broadcast nss"
 915       call flush(iout)
 916       call MPI_Bcast(ihpb_MD(1,1,iprot),maxres*icount,MPI_INTEGER,
 917      &    Master,WHAM_COMM,IERROR)
 918       call MPI_Bcast(jhpb_MD(1,1,iprot),maxres*icount,MPI_INTEGER,
 919      &    Master,WHAM_COMM,IERROR)
 920       write (iout,*) "Finished broadcast ihpb jhpb"
 921       call flush(iout)
 922       call MPI_Bcast(cg_MD(1,1,1,iprot),3*icount*maxres2,
 923      &    MPI_REAL,Master,WHAM_COMM,IERROR)
 924       write (iout,*) "Finished broadcast cg_MD"
 925       call flush(iout)
 926       call MPI_Bcast(cgforce_MD(1,1,1,iprot),3*icount*maxres2,
 927      &    MPI_REAL,Master,WHAM_COMM,IERROR)
 928       write (iout,*) "Finished broadcast cgforce_MD"
 929       call flush(iout)
 930       endif
 931       write (iout,*) "Processor",me," finished broadcast"
 932       call flush(iout)
 933       call dawrite_MDcoords(iprot,jj_old,jj,ientout)
 934 #ifdef DEBUG
 935       write (iout,*) "Protein",iprot
 936       write(iout,*)"Processor",me," received",icount,
 937      &  " MD conformations and forces"
 938       do i=1,icount
 939         write (iout,'(8f10.4)') ((cg_MD(l,k,i,iprot),l=1,3),k=1,nres)
 940         write(iout,'(8f10.4)')((cg_MD(l,k+nres,i,iprot),l=1,3),
 941      &    k=nnt,nct)
 942         write (iout,'(16i5)') nss_MD(i,iprot),(ihpb_MD(k,i,iprot),
 943      &    jhpb_MD(k,i,iprot),k=1,nss_MD(i,iprot))
 944         write(iout,'(8f10.4)')((cgforce_MD(l,k,i,iprot),l=1,3),k=1,
 945      &    nsite(iprot))
 946       enddo
 947 #endif
 948       nconf_forc(iprot)=jj
 949       write (iout,*) "jj",jj," nconf_forc",nconf_forc(iprot)
 950       enddo
 951       return
 952       end
 953 c-------------------------------------------------------------------------------
 954       subroutine work_partition_MD(lprint)
 955 c Split the conformations between processors
 956       implicit none
 957       include "DIMENSIONS"
 958       include "DIMENSIONS.ZSCOPT"
 959       include "mpif.h"
 960       include "COMMON.CLASSES"
 961       include "COMMON.IOUNITS"
 962       include "COMMON.MPI"
 963       include "COMMON.VMCPAR"
 964       include "COMMON.CONTROL"
 965       integer n,chunk,iprot,i,j,ii,remainder
 966       integer kolor,key,ierror,errcode
 967       logical lprint
 968 C
 969 C Divide conformations between processors; for each proteins the first and
 970 C the last conformation to handle by ith processor is stored in
 971 C indstart(i,iprot) and indend(i,iprot), respectively.
 972 C
 973 C First try to assign equal number of conformations to each processor.
 974 C
 975       do iprot=1,nprot
 976
 977         if (.not.fmatch(iprot)) cycle
 978
 979         n=ntot_work_MD(iprot)
 980         write (iout,*) "Protein",iprot," n=",n
 981         indstart_MD(0,iprot)=1
 982         chunk = N/nprocs1
 983         scount_MD(0,iprot) = chunk
 984 c        print *,"i",0," indstart",indstart(0,iprot)," scount",
 985 c     &     scount(0,iprot)
 986         do i=1,nprocs1-1
 987           indstart_MD(i,iprot)=chunk+indstart_MD(i-1,iprot)
 988           scount_MD(i,iprot)=scount_MD(i-1,iprot)
 989 c          print *,"i",i," indstart",indstart(i,iprot)," scount",
 990 c     &     scount(i,iprot)
 991         enddo
 992 C
 993 C Determine how many conformations remained yet unassigned.
 994 C
 995         remainder=N-(indstart_MD(nprocs1-1,iprot)
 996      &    +scount_MD(nprocs1-1,iprot)-1)
 997 c        print *,"remainder",remainder
 998 C
 999 C Assign the remainder conformations to consecutive processors, starting
1000 C from the lowest rank; this continues until the list is exhausted.
1001 C
1002         if (remainder .gt. 0) then
1003           do i=1,remainder
1004             scount_MD(i-1,iprot) = scount_MD(i-1,iprot) + 1
1005             indstart_MD(i,iprot) = indstart_MD(i,iprot) + i
1006           enddo
1007           do i=remainder+1,nprocs1-1
1008             indstart_MD(i,iprot) = indstart_MD(i,iprot) + remainder
1009           enddo
1010         endif
1011
1012         indstart_MD(nprocs1,iprot)=N+1
1013         scount_MD(nprocs1,iprot)=0
1014
1015         do i=0,NProcs1
1016           indend_MD(i,iprot)=indstart_MD(i,iprot)+scount_MD(i,iprot)-1
1017           idispl_MD(i,iprot)=indstart_MD(i,iprot)-1
1018         enddo
1019
1020         N=0
1021         do i=0,Nprocs1-1
1022           N=N+indend_MD(i,iprot)-indstart_MD(i,iprot)+1
1023         enddo
1024
1025 c        print *,"N",n," NTOT",ntot_work(iprot)
1026         if (N.ne.ntot_work_MD(iprot)) then
1027           write (iout,*) "!!! Checksum error on processor",me
1028           call flush(iout)
1029           call MPI_Abort( WHAM_COMM, Ierror, Errcode )
1030         endif
1031
1032       enddo ! iprot
1033       if (lprint) then
1034         write (iout,*) "Partition of work between processors"
1035         do iprot=1,nprot
1036           write (iout,*) "Protein",iprot
1037           do i=0,nprocs1-1
1038             write (iout,'(a,i5,a,i7,a,i7,a,i7)')
1039      &        "Processor",i," indstart_MD",indstart_MD(i,iprot),
1040      &        " indend_MD",indend_MD(i,iprot),
1041      &        " count_MD",scount_MD(i,iprot)
1042           enddo
1043         enddo
1044       endif
1045       return
1046       end
1047 c------------------------------------------------------------------------------
1048 #endif