subroutine atom_partition(iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.IOUNITS" include "COMMON.FMATCH" include "COMMON.INTERACT" include "COMMON.NAMES" include "COMMON.CHAIN" include "COMMON.PROTFILES" character*4 recid,atnam,resnam integer iprot,ires_prev,ires,nnres,nsc,iat,i,ii,j integer ilen external ilen logical carflag write (iout,*) "In atom_partition" call restore_molinfo(iprot) write (iout,*) "nres",nres open(ipdbin,file=fpdbfile(iprot),status="old",err=10) goto 11 10 write (iout,*) "Error opening all-atom-type file ", & fpdbfile(iprot)(:ilen(fpdbfile(iprot))) stop 11 continue write (iout,*) "atom_partitio: after opening file" #ifdef PISIORNIA do i=1,nres-1 print *,"i",i read(ipdbin,'(i3,5x,i3,1x,$)',err=20,end=20) iii,natp(i,iprot) do j=1,natp(i,iprot) read(ipdbin,'(a4,i4,1x,$)',err=20,end=20) atp(j,i,iprot), & iatp(j,i,iprot) enddo read(ipdbin,'(3x,i3,1x,$)',err=20,end=20) natsc(i,iprot) do j=1,natsc(i,iprot) read(ipdbin,'(a4,i4,1x,$)',err=20,end=20) atsc(j,i,iprot), & iatsc(j,i,iprot) enddo read(ipdbin,*) enddo #else natp(:,iprot)=0 natsc(:,iprot)=1 natsc(1,iprot)=0 natsc(nres,iprot)=0 ires_prev=0 nat(iprot)=0 nsc=0 do read (ipdbin,'(a4,i7,1x,a4,1x,a3,i6)',end=20,err=20) & recid,iat,atnam,resnam,ires write (iout,*) recid,iat,atnam,resnam,ires if (recid.ne."ATOM") cycle if (iat.gt.nat(iprot)) nat(iprot)=iat if (ires.ne.ires_prev) then carflag=.false. ires_prev=ires endif if (atnam.eq." CA ") then if (restyp(itype(ires+nnt-1)).ne.resnam) then write (iout,*) & "Error: different residue name in force-match template", & " protein",iprot," residue",ires,resnam,restyp(itype(ires+1)) stop endif carflag=.true. iatsc(1,ires+nnt-1,iprot)=iat atsc(1,ires+nnt-1,iprot)=atnam(2:) endif c print *,"carflag ",carflag if (.not.carflag) then if (resnam.eq."PRO ") then write (iout,*) resnam,ires c if (atnam.eq." N ".or.atnam.eq." CD ".or.atnam.eq." HD2" c if (atnam.eq." N ".or.atnam.eq." HD2".or.atnam.eq." HD3") if (atnam.eq." N ") & then natp(ires-1+nnt-1,iprot)=natp(ires-1+nnt-1,iprot)+1 c print *,"back",natp(ires-1) ii=natp(ires-1+nnt-1,iprot) atp(ii,ires-1+nnt-1,iprot)=atnam(2:) iatp(ii,ires-1+nnt-1,iprot)=iat else natsc(ires+nnt-1,iprot)=natsc(ires+nnt-1,iprot)+1 c print *,"sc",natsc(ires) ii=natsc(ires+nnt-1,iprot) iatsc(ii,ires+nnt-1,iprot)=iat atsc(ii,ires+nnt-1,iprot)=atnam(2:) endif else natp(ires-1+nnt-1,iprot)=natp(ires-1+nnt-1,iprot)+1 ii=natp(ires-1+nnt-1,iprot) atp(ii,ires-1+nnt-1,iprot)=atnam(2:) iatp(ii,ires-1+nnt-1,iprot)=iat endif elseif(atnam.eq." C ".or.atnam.eq." O " .or.atnam.eq." OXT") & then natp(ires+nnt-1,iprot)=natp(ires+nnt-1,iprot)+1 ii=natp(ires+nnt-1,iprot) if (atnam(1:1).eq." ") then atp(ii,ires+nnt-1,iprot)=atnam(2:) else atp(ii,ires+nnt-1,iprot)=atnam endif iatp(ii,ires+nnt-1,iprot)=iat else if (atnam.ne. " CA ") then natsc(ires+nnt-1,iprot)=natsc(ires+nnt-1,iprot)+1 ii=natsc(ires+nnt-1,iprot) iatsc(ii,ires+nnt-1,iprot)=iat if (atnam(1:1).eq." ") then atsc(ii,ires+nnt-1,iprot)=atnam(2:) else atsc(ii,ires+nnt-1,iprot)=atnam endif endif nnres=ires if (atnam.eq." CA ".and.resnam.ne."GLY ") nsc=nsc+1 enddo 20 continue close(ipdbin) write (iout,*) "End reading nnres",nnres," nres",nres, & " nsc",nsc," nat",nat(iprot) #endif close(ipdbin) c--------------------------------------------------------------------- do i=1,nres write(iout,'(i3,1x,a3,1x,5h back,i3,1x,$)')i,restyp(itype((i))), & natp(i,iprot) do j=1,natp(i,iprot) write(iout,'(a4,i4,1x,$)') atp(j,i,iprot),iatp(j,i,iprot) enddo write(iout,'(i3,1x,3h sc,i3,1x,$)') i,natsc(i,iprot) do j=1,natsc(i,iprot) write(iout,'(a4,i4,1x,$)') atsc(j,i,iprot),iatsc(j,i,iprot) enddo write (iout,*) enddo nsite(iprot)=nres do i=nnt,nct if (itype(i).ne.10) nsite(iprot)=nsite(iprot)+1 enddo write (iout,*) "iprot",iprot," nsite",nsite(iprot),nres+nsc return end c--------------------------------------------------------------------- subroutine read_allat_coord_forces(iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror,errcode #endif include "COMMON.CONTROL" include "COMMON.IOUNITS" include "COMMON.FMATCH" include "COMMON.CHAIN" include "COMMON.INTERACT" include "COMMON.PROTFILES" include "COMMON.SBRIDGE" include "COMMON.NAMES" include "COMMON.WEIGHTS" include "COMMON.CLASSES" include "COMMON.PROTNAME" integer i,ii,j,k,l,kkk,kkkk,jj_old,jj,jjj,icount,iis,nnsc integer iprot integer ilen external ilen double precision casc(3),dss jj_old=1 call restore_molinfo(iprot) open (ientout,file=bprotfiles_MD(iprot),status="unknown", & form="unformatted",access="direct",recl=lenrec_MD(iprot)) c Change AL 12/30/2017 #ifdef MPI if (me.eq.Master) then ! Loop over snapshots #endif! open(ipdbin,file=fcoordfile(iprot),status="old",err=10) goto 11 10 write (iout,*) "Error opening MD coordinate file ", & fcoordfile(iprot)(:ilen(fcoordfile(iprot))) call flush(iout) #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR) #else stop #endif 11 continue open(ientin,file=fforcefile(iprot),status="old",err=20) goto 21 20 write (iout,*) "Error opening MD force file ", & fforcefile(iprot)(:ilen(fforcefile(iprot))) call flush(iout) #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR) #else stop #endif 21 continue write (iout,*) "nat",nat(iprot) call flush(iout) kkkk=0 jjj=0 jj=0 icount=0 read (ipdbin,*) read (ientin,*) do kkk=1,maxstr c do kkk=1,1 read (ipdbin,'(10f8.3)',end=30,err=30) & ((atcoord(j,i),j=1,3),i=1,nat(iprot)) read(ientin,'(10f8.3)',end=30,err=30)((atforce(j,i),j=1,3),i=1, & nat(iprot)) if (mod(kkk,ifreq_MD(iprot)).gt.0) cycle kkkk = kkkk + 1 #ifdef DEBUG write(iout,*) "All-atom Coordinates and forces",kkk,kkkk do i=1,nres-1 do j=1,natp(i,iprot) write(iout,'(2i3,1x,a4,i5,3f8.3,5x,3f8.3)') & i,j,atp(j,i,iprot),iatp(j,i,iprot), & (atcoord(k,iatp(j,i,iprot)),k=1,3), & (atforce(k,iatp(j,i,iprot)),k=1,3) enddo do j=1,natsc(i,iprot) write(iout,'(2i3,1x,a4,i5,3f8.3,5x,3f8.3)') & i,j,atsc(j,i,iprot),iatsc(j,i,iprot), & (atcoord(k,iatsc(j,i,iprot)),k=1,3), & (atforce(k,iatsc(j,i,iprot)),k=1,3) enddo enddo #endif do k=1,3 c(k,1)=atcoord(k,iatp(1,1,iprot)) c(k,nres)=atcoord(k,iatp(3,nres-1,iprot)) enddo do i=nnt,nct do j=1,3 c(j,i)=atcoord(j,iatsc(1,i,iprot)) enddo enddo iis=0 do i=1,nres casc=0.0d0 c write (iout,*) "i",i," natsc",natsc(i) nnsc=0 do j=1,natsc(i,iprot) c write (iout,*) "j",j," iatsc",iatsc(j,i) if (atsc(j,i,iprot)(1:1).eq."H") cycle nnsc=nnsc+1 do k=1,3 casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot)) enddo enddo do k=1,3 c(k,i+nres)=casc(k)/nnsc enddo c write (iout,*) (c(k,i+nres),k=1,3) enddo c Detect disulfide bonds nss=0 do i=nnt,nct do j=i+1,nct if (itype(i).eq.1 .and. itype(j).eq.1) then do k=1,natsc(i,iprot) do l=1,natsc(j,iprot) if (atsc(k,i,iprot).eq." S ".and. & atsc(l,j,iprot).eq." S ") then ii=iatsc(k,i,iprot) jj=iatsc(l,j,iprot) dSS=dsqrt((atcoord(1,ii)-atcoord(1,jj))**2+ & (atcoord(2,ii)-atcoord(2,jj))**2+ & (atcoord(3,ii)-atcoord(3,jj))**2) endif if (dSS.lt.2.5d0) then nss=nss+1 ihpb(nss)=ii jhpb(nss)=jj endif enddo enddo endif enddo enddo #ifdef DEBUG write (iout,*) "CG Coordinates" do i=1,nres if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then write (iout,'(a4,i5,3f10.3)')restyp(itype(i)),i,(c(j,i),j=1,3) else write (iout,'(a4,i5,3f10.3,5x,3f10.3)') & restyp(itype(i)),i,(c(j,i),j=1,3),(c(j,i+nres),j=1,3) endif enddo write (iout,'(8f10.5)') & ((c(l,k),l=1,3),k=1,nres), & ((c(l,k+nres),l=1,3),k=nnt,nct) call int_from_cart1(.true.) #endif call add_new_MDconf(jjj,jj,jj_old,icount,iprot) enddo 30 continue call write_and_send_MDconf(icount,jj_old,jj,iprot) nconf_forc(iprot)=kkkk ntot_work_MD(iprot)=kkkk call write_and_send_MDconf(0,jj_old,jj,iprot) #ifdef MPI else #endif jjj=0 jj=0 icount=0 do i=1,maxstr list_conf_MD(i,iprot)=i enddo call receive_and_write_MDconf(icount,jj_old,jj,iprot) #ifdef MPI endif close(icbase) ntot_MD(iprot)=jj write (iout,*) "Protein", & protname(iprot)(:ilen(protname(iprot))),", ",ntot_MD(iprot), & " conformatons read ",ntot_MD(iprot), & " conformations written to ", & bprotfiles_MD(iprot)(:ilen(bprotfiles_MD(iprot))) ntot_MD(0) = ntot_MD(0)+ntot_MD(iprot) #endif return end c------------------------------------------------------------------- subroutine compute_CG_forces(kkkk,iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.IOUNITS" include "COMMON.FMATCH" include "COMMON.CHAIN" include "COMMON.INTERACT" include "COMMON.NAMES" integer i,ii,iis,j,k,kkkk,iprot,nsc,nscc,nnsc,n,nn double precision yy(3),caca(3),casc(3),xx,cacanorm,cascnorm double precision scalar ii=0 do j=1,natp(1,iprot) c write (iout,'(3i5,3f10.5)') kkkk,j,iatp(j,1,iprot), c & (atforce(k,iatp(j,1,iprot)),k=1,3) do k=1,3 CGforce_MD(k,nnt,kkkk,iprot)=CGforce_MD(k,nnt,kkkk,iprot) & +atforce(k,iatp(j,1,iprot)) enddo enddo do i=nnt,nct-1 do j=1,3 caca(j)=atcoord(j,iatsc(1,i+1,iprot)) & -atcoord(j,iatsc(1,i,iprot)) enddo cacanorm=dsqrt(scalar(caca,caca)) c write(iout,*)"i",i," cacanorm",cacanorm do k=1,3 caca(k)=caca(k)/cacanorm enddo do j=1,natp(i,iprot) yy=0.0d0 do k=1,3 yy(k)=atcoord(k,iatp(j,i,iprot))-atcoord(k,iatsc(1,i,iprot)) enddo xx=scalar(caca,yy)/cacanorm do k=1,3 CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot) & +(1-xx)*atforce(k,iatp(j,i,iprot)) CGforce_MD(k,i+1,kkkk,iprot)=CGforce_MD(k,i+1,kkkk,iprot) & +xx*atforce(k,iatp(j,i,iprot)) enddo enddo enddo c#ifdef SC nscc=0 iis=nres do i=nnt,nct casc=0.0d0 c print *,"i",i," natsc",natsc(i) nnsc=0 do j=1,natsc(i,iprot) c print *,"j",j," iatsc",iatsc(j,i) if (atsc(j,i,iprot)(1:1).eq."H") cycle nnsc=nnsc+1 do k=1,3 casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot)) enddo enddo do k=1,3 casc(k)=casc(k)/nnsc-atcoord(k,iatsc(1,i,iprot)) enddo cascnorm=dsqrt(scalar(casc,casc)) c write(iout,*)"i",i,restyp(itype(i))," natsc",natsc(i,iprot), c & " nnsc",nnsc," cascnorm",cascnorm if (nnsc.gt.1) then nscc=nscc+1 iis=iis+1 ires_sc(iis,iprot)=i endif c write (iout,*) "nnsc",nnsc," iis",iis if (cascnorm.gt.1.0d-2) then do k=1,3 casc(k)=casc(k)/cascnorm enddo else casc=0.0d0 endif do k=1,3 CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot) & +atforce(k,iatsc(1,i,iprot)) enddo do j=2,natsc(i,iprot) if (atsc(j,i,iprot)(1:2).eq."HA") then do k=1,3 CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot) & +atforce(k,iatsc(j,i,iprot)) enddo c write (4,*) "HA",i,(CGforce_MD(k,i,kkkk),k=1,3) else do k=1,3 yy(k)=atcoord(k,iatsc(j,i,iprot))-atcoord(k,iatsc(1,i,iprot)) enddo xx=scalar(casc,yy)/cascnorm do k=1,3 CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot) & +(1.0d0-xx)*atforce(k,iatsc(j,i,iprot)) CGforce_MD(k,iis,kkkk,iprot)=CGforce_MD(k,iis,kkkk,iprot) & +xx*atforce(k,iatsc(j,i,iprot)) enddo c write (4,*)"REG CA",i,(CGforce_MD(k,i,kkkk),k=1,3) c write (4,*)"REG SC",iis,(CGforce_MD(k,iis,kkkk),k=1,3) endif enddo enddo do j=1,natp(nct,iprot) do k=1,3 CGforce_MD(k,nct,kkkk,iprot)=CGforce_MD(k,nct,kkkk,iprot) & +atforce(k,iatp(j,nct,iprot)) enddo enddo c#endif nn = ii c print *,"nres",nres," nsc",nsc n = nres+nsc #ifdef DEBUG write (iout,*) "CG forces calculated from all-atom forces" write (iout,"(4x,a,t55,a)") "Backbone","Sidechain" ii=nres do i=1,nres if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then write (iout,'(a4,i5,3e15.5)') restyp(itype(i)),i, & (CGforce_MD(j,i,kkkk,iprot),j=1,3) else ii=ii+1 write (iout,'(a4,i5,3e15.5,5x,3e15.5)') restyp(itype(i)),i, & (CGforce_MD(j,i,kkkk,iprot),j=1,3), & (CGforce_MD(j,ii,kkkk,iprot),j=1,3) endif enddo #endif return end c------------------------------------------------------------------- subroutine compute_CG_forces_ave(kkkk,iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.IOUNITS" include "COMMON.FMATCH" include "COMMON.CHAIN" include "COMMON.INTERACT" include "COMMON.NAMES" integer i,ii,iis,j,k,kkkk,iprot,nsc,nscc,nnsc,n,nn double precision yy(3),caca(3),casc(3),xx,cacanorm,cascnorm double precision scalar ii=0 do j=1,natp(1,iprot) c write (iout,'(3i5,3f10.5)') kkkk,j,iatp(j,1,iprot), c & (atforce(k,iatp(j,1,iprot)),k=1,3) do k=1,3 CGforce_MD_ave(k,nnt,kkkk,iprot)=CGforce_MD(k,nnt,kkkk,iprot) & +atforce(k,iatp(j,1,iprot)) enddo enddo do i=nnt,nct-1 do j=1,3 caca(j)=atcoord(j,iatsc(1,i+1,iprot)) & -atcoord(j,iatsc(1,i,iprot)) enddo cacanorm=dsqrt(scalar(caca,caca)) c write(iout,*)"i",i," cacanorm",cacanorm do k=1,3 caca(k)=caca(k)/cacanorm enddo do j=1,natp(i,iprot) yy=0.0d0 do k=1,3 yy(k)=atcoord(k,iatp(j,i,iprot))-atcoord(k,iatsc(1,i,iprot)) enddo xx=scalar(caca,yy)/cacanorm do k=1,3 CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot) & +(1-xx)*atforce(k,iatp(j,i,iprot)) CGforce_MD_ave(k,i+1,kkkk,iprot)=CGforce_MD(k,i+1,kkkk,iprot) & +xx*atforce(k,iatp(j,i,iprot)) enddo enddo enddo c#ifdef SC nscc=0 iis=nres do i=nnt,nct casc=0.0d0 c print *,"i",i," natsc",natsc(i) nnsc=0 do j=1,natsc(i,iprot) c print *,"j",j," iatsc",iatsc(j,i) if (atsc(j,i,iprot)(1:1).eq."H") cycle nnsc=nnsc+1 do k=1,3 casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot)) enddo enddo do k=1,3 casc(k)=casc(k)/nnsc-atcoord(k,iatsc(1,i,iprot)) enddo cascnorm=dsqrt(scalar(casc,casc)) c write(iout,*)"i",i,restyp(itype(i))," natsc",natsc(i,iprot), c & " nnsc",nnsc," cascnorm",cascnorm if (nnsc.gt.1) then nscc=nscc+1 iis=iis+1 ires_sc(iis,iprot)=i endif c write (iout,*) "nnsc",nnsc," iis",iis if (cascnorm.gt.1.0d-2) then do k=1,3 casc(k)=casc(k)/cascnorm enddo else casc=0.0d0 endif do k=1,3 CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot) & +atforce(k,iatsc(1,i,iprot)) enddo do j=2,natsc(i,iprot) if (atsc(j,i,iprot)(1:2).eq."HA") then do k=1,3 CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot) & +atforce(k,iatsc(j,i,iprot)) enddo c write (4,*) "HA",i,(CGforce_MD(k,i,kkkk),k=1,3) else do k=1,3 yy(k)=atcoord(k,iatsc(j,i,iprot))-atcoord(k,iatsc(1,i,iprot)) enddo xx=scalar(casc,yy)/cascnorm do k=1,3 CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot) & +(1.0d0-xx)*atforce(k,iatsc(j,i,iprot)) CGforce_MD_ave(k,iis,kkkk,iprot)=CGforce_MD(k,iis,kkkk,iprot) & +xx*atforce(k,iatsc(j,i,iprot)) enddo c write (4,*)"REG CA",i,(CGforce_MD(k,i,kkkk),k=1,3) c write (4,*)"REG SC",iis,(CGforce_MD(k,iis,kkkk),k=1,3) endif enddo enddo do j=1,natp(nct,iprot) do k=1,3 CGforce_MD_ave(k,nct,kkkk,iprot)=CGforce_MD(k,nct,kkkk,iprot) & +atforce(k,iatp(j,nct,iprot)) enddo enddo c#endif nn = ii c print *,"nres",nres," nsc",nsc n = nres+nsc #ifdef DEBUG write (iout,*) "CG forces calculated from all-atom forces" write (iout,"(4x,a,t55,a)") "Backbone","Sidechain" ii=nres do i=1,nres if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then write (iout,'(a4,i5,3e15.5)') restyp(itype(i)),i, & (CGforce_MD_ave(j,i,kkkk,iprot),j=1,3) else ii=ii+1 write (iout,'(a4,i5,3e15.5,5x,3e15.5)') restyp(itype(i)),i, & (CGforce_MD_ave(j,i,kkkk,iprot),j=1,3), & (CGforce_MD_ave(j,ii,kkkk,iprot),j=1,3) endif enddo #endif return end c---------------------------------------------------------------------- subroutine add_new_MDconf(jjj,jj,jj_old,icount,iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CHAIN" include "COMMON.INTERACT" include "COMMON.LOCAL" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.PROTNAME" include "COMMON.VMCPAR" include "COMMON.WEIGHTS" include "COMMON.NAMES" include "COMMON.ALLPROT" include "COMMON.WEIGHTDER" include "COMMON.VAR" include "COMMON.SBRIDGE" include "COMMON.GEO" include "COMMON.COMPAR" integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,ib,ibatch, & nn,nn1,inan,Next,itj jjj=jjj+1 c write(iout,*) "add_new_MDconf jjj jj jj_old:",jjj,jj,jj_old call flush(iout) call int_from_cart1(.false.) call flush(iout) c write (iout,*) "After int_from_cart1" call flush(iout) do j=nnt+1,nct if (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.5d0) then write (iout,*) "nnt",nnt," nct",nct write (iout,*) "Bad CA-CA bond length",j," ",vbld(j) write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j) write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) & "This conformation WILL NOT be added to the database." return endif enddo do j=nnt,nct itj=itype(j) if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then write (iout,*) "nnt",nnt," nct",nct write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j) write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) & "This conformation WILL NOT be added to the database." return endif enddo do j=3,nres if (theta(j).le.0.0d0) then write (iout,*) & "Zero theta angle(s) in conformation",ii write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) & "This conformation WILL NOT be added to the database." return endif if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad enddo jj=jj+1 icount=icount+1 list_conf_MD(jj,iprot)=jj c write (iout,*) "calling store_MDconf",icount call flush(iout) call store_MDconf_from_file(jj,icount,iprot) c write (iout,*) "finished store_MDconf" call flush(iout) if (icount.eq.maxstr_proc) then #ifdef DEBUG write (iout,* ) "jj_old",jj_old," jj",jj write (iout,*) "Calling write_and_send_cconf" call flush(iout) #endif call write_and_send_MDconf(icount,jj_old,jj,iprot) jj_old=jj+1 #ifdef DEBUG write (iout,*) "Exited write_and_send_cconf" call flush(iout) #endif icount=0 endif return end c------------------------------------------------------------------------------ subroutine store_MDconf_from_file(jj,icount,iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CHAIN" include "COMMON.SBRIDGE" include "COMMON.INTERACT" include "COMMON.IOUNITS" include "COMMON.CLASSES" include "COMMON.ALLPROT" include "COMMON.VAR" include "COMMON.FMATCH" integer i,j,jj,icount,ibatch,iprot c Store the conformation that has been read in do i=1,2*nres do j=1,3 cg_MD(j,i,icount,iprot)=c(j,i) enddo enddo nss_MD(icount,iprot)=nss do i=1,nss ihpb_MD(i,icount,iprot)=ihpb(i) jhpb_MD(i,icount,iprot)=jhpb(i) enddo call compute_CG_forces(icount,iprot) return end c------------------------------------------------------------------------------ subroutine write_and_send_MDconf(icount,jj_old,jj,iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR include "COMMON.MPI" #endif include "COMMON.WEIGHTS" include "COMMON.CHAIN" include "COMMON.SBRIDGE" include "COMMON.INTERACT" include "COMMON.IOUNITS" include "COMMON.CLASSES" include "COMMON.VAR" include "COMMON.ALLPROT" include "COMMON.ENERGIES" include "COMMON.WEIGHTDER" include "COMMON.OPTIM" include "COMMON.COMPAR" include "COMMON.FMATCH" integer icount,jj_old,jj,Next,iprot c Write the structures to a scratch file #ifdef MPI #ifdef DEBUG write (iout,*) "Processor",me," entered WRITE_AND_SEND_MDCONF" write (iout,*) "iprot",iprot call flush(iout) #endif call MPI_Bcast(icount,1,MPI_INTEGER,Master,WHAM_COMM,IERROR) #ifdef DEBUG write (iout,*) "Processor",me," sent icount=",icount write (iout,*) "Processor",me," jj_old",jj_old," jj",jj call flush(iout) #endif if (icount.gt.0) then call MPI_Bcast(jj_old,1,MPI_INTEGER,Master,WHAM_COMM,IERROR) call MPI_Bcast(jj,1,MPI_INTEGER,Master,WHAM_COMM,IERROR) write (iout,*) "Finished broadcast jj_old jj" call flush(iout) call MPI_Bcast(nss_MD(1,iprot),icount,MPI_INTEGER,Master, & WHAM_COMM,IERROR) write (iout,*) "Finished broadcast nss" call flush(iout) call MPI_Bcast(ihpb_MD(1,1,iprot),maxres*icount,MPI_INTEGER, & Master,WHAM_COMM,IERROR) call MPI_Bcast(jhpb_MD(1,1,iprot),maxres*icount,MPI_INTEGER, & Master,WHAM_COMM,IERROR) write (iout,*) "Finished broadcast ihpb jhpb" call flush(iout) call MPI_Bcast(cg_MD(1,1,1,iprot),3*icount*maxres2, & MPI_REAL,Master,WHAM_COMM,IERROR) write (iout,*) "Finished broadcast cg_MD" call flush(iout) call MPI_Bcast(cgforce_MD(1,1,1,iprot),3*icount*maxres2, & MPI_REAL,Master,WHAM_COMM,IERROR) write (iout,*) "Finished broadcast cgforce_MD" call flush(iout) endif write (iout,*) "Master finished broadcast" call flush(iout) #endif call dawrite_MDcoords(iprot,jj_old,jj,ientout) nconf_forc(iprot)=jj write (iout,*) "iprot",iprot," jj",jj," nconf_forc", & nconf_forc(iprot) #ifdef DEBUG write (iout,*) "Protein",iprot write(iout,*)"Processor",me," received",icount, & " MD conformations and forces" do i=1,icount write (iout,'(8f10.4)') ((cg_MD(l,k,i,iprot),l=1,3,k=1,nres) write(iout,'(8f10.4)')((cg_MD(l,k,i+nres,iprot),l=1,3,k=nnt,nct) write (iout,'(16i5)') nss_MD(i,iprot),(ihpb_MD(k,i,iprot), & jhpb_MD(k,i,iprot),k=1,nss_MD(i,iprot)) write(iout,'(8f10.4)')((cgforce_MD(l,k,i,iprot),l=1,3,k=,1, & nsite(iprot))) enddo #endif return end c------------------------------------------------------------------------------ #ifdef MPI subroutine receive_and_write_MDconf(icount,jj_old,jj,iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "mpif.h" integer IERROR,STATUS(MPI_STATUS_SIZE) include "COMMON.MPI" include "COMMON.CHAIN" include "COMMON.SBRIDGE" include "COMMON.INTERACT" include "COMMON.IOUNITS" include "COMMON.CLASSES" include "COMMON.FMATCH" include "COMMON.ENERGIES" include "COMMON.VAR" include "COMMON.GEO" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.COMPAR" include "COMMON.OPTIM" integer i,j,k,icount,jj_old,jj,iprot,Previous,Next icount=1 do while (icount.gt.0) #ifdef DEBUG write (iout,*) "Processor",me," entered RECEIVE_AND_SEND_MDCONF" write (iout,*) "iprot",iprot call flush(iout) #endif call MPI_Bcast(icount,1,MPI_INTEGER,Master,WHAM_COMM,IERROR) #ifdef DEBUG write (iout,*) "Processor",me," received icount=",icount call flush(iout) #endif if (icount.gt.0) then call MPI_Bcast(jj_old,1,MPI_INTEGER,Master,WHAM_COMM,IERROR) call MPI_Bcast(jj,1,MPI_INTEGER,Master,WHAM_COMM,IERROR) write (iout,*) "Finished broadcast jj_old jj" write (iout,*) "Processor",me," jj_old",jj_old," jj",jj call flush(iout) call MPI_Bcast(nss_MD(1,iprot),icount,MPI_INTEGER,Master, & WHAM_COMM,IERROR) write (iout,*) "Finished broadcast nss" call flush(iout) call MPI_Bcast(ihpb_MD(1,1,iprot),maxres*icount,MPI_INTEGER, & Master,WHAM_COMM,IERROR) call MPI_Bcast(jhpb_MD(1,1,iprot),maxres*icount,MPI_INTEGER, & Master,WHAM_COMM,IERROR) write (iout,*) "Finished broadcast ihpb jhpb" call flush(iout) call MPI_Bcast(cg_MD(1,1,1,iprot),3*icount*maxres2, & MPI_REAL,Master,WHAM_COMM,IERROR) write (iout,*) "Finished broadcast cg_MD" call flush(iout) call MPI_Bcast(cgforce_MD(1,1,1,iprot),3*icount*maxres2, & MPI_REAL,Master,WHAM_COMM,IERROR) write (iout,*) "Finished broadcast cgforce_MD" call flush(iout) endif write (iout,*) "Processor",me," finished broadcast" call flush(iout) call dawrite_MDcoords(iprot,jj_old,jj,ientout) #ifdef DEBUG write (iout,*) "Protein",iprot write(iout,*)"Processor",me," received",icount, & " MD conformations and forces" do i=1,icount write (iout,'(8f10.4)') ((cg_MD(l,k,i,iprot),l=1,3,k=1,nres) write(iout,'(8f10.4)')((cg_MD(l,k,i+nres,iprot),l=1,3,k=nnt,nct) write (iout,'(16i5)') nss_MD(i,iprot),(ihpb_MD(k,i,iprot), & jhpb_MD(k,i,iprot),k=1,nss_MD(i,iprot)) write(iout,'(8f10.4)')((cgforce_MD(l,k,i,iprot),l=1,3,k=,1, & nsite(iprot))) enddo #endif nconf_forc(iprot)=jj write (iout,*) "jj",jj," nconf_forc",nconf_forc(iprot) enddo return end c------------------------------------------------------------------------------- subroutine work_partition_MD(lprint) c Split the conformations between processors implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "mpif.h" include "COMMON.CLASSES" include "COMMON.IOUNITS" include "COMMON.MPI" include "COMMON.VMCPAR" include "COMMON.CONTROL" integer n,chunk,iprot,i,j,ii,remainder integer kolor,key,ierror,errcode logical lprint C C Divide conformations between processors; for each proteins the first and C the last conformation to handle by ith processor is stored in C indstart(i,iprot) and indend(i,iprot), respectively. C C First try to assign equal number of conformations to each processor. C do iprot=1,nprot if (.not.fmatch(iprot)) cycle n=ntot_work_MD(iprot) write (iout,*) "Protein",iprot," n=",n indstart_MD(0,iprot)=1 chunk = N/nprocs1 scount_MD(0,iprot) = chunk c print *,"i",0," indstart",indstart(0,iprot)," scount", c & scount(0,iprot) do i=1,nprocs1-1 indstart_MD(i,iprot)=chunk+indstart_MD(i-1,iprot) scount_MD(i,iprot)=scount_MD(i-1,iprot) c print *,"i",i," indstart",indstart(i,iprot)," scount", c & scount(i,iprot) enddo C C Determine how many conformations remained yet unassigned. C remainder=N-(indstart_MD(nprocs1-1,iprot) & +scount_MD(nprocs1-1,iprot)-1) c print *,"remainder",remainder C C Assign the remainder conformations to consecutive processors, starting C from the lowest rank; this continues until the list is exhausted. C if (remainder .gt. 0) then do i=1,remainder scount_MD(i-1,iprot) = scount_MD(i-1,iprot) + 1 indstart_MD(i,iprot) = indstart_MD(i,iprot) + i enddo do i=remainder+1,nprocs1-1 indstart_MD(i,iprot) = indstart_MD(i,iprot) + remainder enddo endif indstart_MD(nprocs1,iprot)=N+1 scount_MD(nprocs1,iprot)=0 do i=0,NProcs1 indend_MD(i,iprot)=indstart_MD(i,iprot)+scount_MD(i,iprot)-1 idispl_MD(i,iprot)=indstart_MD(i,iprot)-1 enddo N=0 do i=0,Nprocs1-1 N=N+indend_MD(i,iprot)-indstart_MD(i,iprot)+1 enddo c print *,"N",n," NTOT",ntot_work(iprot) if (N.ne.ntot_MD(iprot)) then write (iout,*) "!!! Checksum error on processor",me call flush(iout) call MPI_Abort( WHAM_COMM, Ierror, Errcode ) endif enddo ! iprot if (lprint) then write (iout,*) "Partition of work between processors" do iprot=1,nprot write (iout,*) "Protein",iprot do i=0,nprocs1-1 write (iout,'(a,i5,a,i7,a,i7,a,i7)') & "Processor",i," indstart_MD",indstart_MD(i,iprot), & " indend_MD",indend_MD(i,iprot), & " count_MD",scount_MD(i,iprot) enddo enddo endif return end c------------------------------------------------------------------------------ #endif