source/maxlik/src-Fmatch_safe/readrtns_forcefit.F.safe

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