1 subroutine atom_partition(iprot)
4 include "DIMENSIONS.ZSCOPT"
5 include "COMMON.IOUNITS"
6 include "COMMON.FMATCH"
7 include "COMMON.INTERACT"
10 include "COMMON.PROTFILES"
11 character*4 recid,atnam,resnam
12 integer iprot,ires_prev,ires,nnres,nsc,iat,i,ii,j
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)
21 10 write (iout,*) "Error opening all-atom-type file ",
22 & fpdbfile(iprot)(:ilen(fpdbfile(iprot)))
25 write (iout,*) "atom_partitio: after opening file"
29 read(ipdbin,'(i3,5x,i3,1x,$)',err=20,end=20) iii,natp(i,iprot)
31 read(ipdbin,'(a4,i4,1x,$)',err=20,end=20) atp(j,i,iprot),
34 read(ipdbin,'(3x,i3,1x,$)',err=20,end=20) natsc(i,iprot)
36 read(ipdbin,'(a4,i4,1x,$)',err=20,end=20) atsc(j,i,iprot),
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
59 if (atnam.eq." CA ") then
60 if (restyp(itype(ires+nnt-1)).ne.resnam) then
62 & "Error: different residue name in force-match template",
63 & " protein",iprot," residue",ires,resnam,restyp(itype(ires+1))
67 iatsc(1,ires+nnt-1,iprot)=iat
68 atsc(1,ires+nnt-1,iprot)=atnam(2:)
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")
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
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:)
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
96 elseif(atnam.eq." C ".or.atnam.eq." O " .or.atnam.eq." OXT")
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:)
103 atp(ii,ires+nnt-1,iprot)=atnam
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:)
113 atsc(ii,ires+nnt-1,iprot)=atnam
117 if (atnam.eq." CA ".and.resnam.ne."GLY ") nsc=nsc+1
121 write (iout,*) "End reading nnres",nnres," nres",nres,
122 & " nsc",nsc," nat",nat(iprot)
125 c---------------------------------------------------------------------
127 write(iout,'(i3,1x,a3,1x,5h back,i3,1x,$)')i,restyp(itype((i))),
130 write(iout,'(a4,i4,1x,$)') atp(j,i,iprot),iatp(j,i,iprot)
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)
140 if (itype(i).ne.10) nsite(iprot)=nsite(iprot)+1
142 write (iout,*) "iprot",iprot," nsite",nsite(iprot),nres+nsc
145 c---------------------------------------------------------------------
146 subroutine read_allat_coord_forces(iprot)
149 include "DIMENSIONS.ZSCOPT"
153 integer ierror,errcode
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
170 double precision casc(3),dss
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
179 if (me.eq.Master) then
180 ! Loop over snapshots
182 open(ipdbin,file=fcoordfile(iprot),status="old",err=10)
184 10 write (iout,*) "Error opening MD coordinate file ",
185 & fcoordfile(iprot)(:ilen(fcoordfile(iprot)))
188 call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR)
193 open(ientin,file=fforcefile(iprot),status="old",err=20)
195 20 write (iout,*) "Error opening MD force file ",
196 & fforcefile(iprot)(:ilen(fforcefile(iprot)))
199 call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR)
204 write (iout,*) "nat",nat(iprot)
210 CGforce_MD_ave(:,:,:,iprot)=0.0d0
211 write (iout,*) "istart_MD",istart_MD(iprot),
212 & " iend_MD",iend_MD(iprot)
215 do kkk=1,iend_MD(iprot)
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,
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
228 CGforce_MD_ave(j,i,kkkk,iprot)=CGforce_MD_ave(j,i,kkkk,iprot)/
233 write(iout,*) "All-atom Coordinates and forces",kkk,kkkk
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)
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)
250 c(k,1)=atcoord(k,iatp(1,1,iprot))
251 c(k,nres)=atcoord(k,iatp(3,nres-1,iprot))
255 c(j,i)=atcoord(j,iatsc(1,i,iprot))
261 c write (iout,*) "i",i," natsc",natsc(i)
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
268 casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot))
272 c(k,i+nres)=casc(k)/nnsc
274 c write (iout,*) (c(k,i+nres),k=1,3)
276 c Detect disulfide bonds
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
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)
291 if (dSS.lt.2.5d0) then
302 write (iout,*) "CG Coordinates"
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)
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)
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.)
317 call add_new_MDconf(jjj,jj,jj_old,icount,iprot)
319 write (iout,*) "Average and momentary CG forces",
320 & " calculated from all-atom forces"
321 write (iout,"(4x,a,t55,a)") "Backbone","Sidechain"
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)
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)
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
348 call write_and_send_MDconf(0,jj_old,jj,iprot)
356 list_conf_MD(i,iprot)=i
358 call receive_and_write_MDconf(icount,jj_old,jj,iprot)
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)
372 c-------------------------------------------------------------------
373 subroutine compute_CG_forces(kkkk,iprot)
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
386 CGforce_MD(:,:,kkkk,iprot) = 0.0d0
388 c write (iout,'(3i5,3f10.5)') kkkk,j,iatp(j,1,iprot),
389 c & (atforce(k,iatp(j,1,iprot)),k=1,3)
391 CGforce_MD(k,nnt,kkkk,iprot)=CGforce_MD(k,nnt,kkkk,iprot)
392 & +atforce(k,iatp(j,1,iprot))
397 caca(j)=atcoord(j,iatsc(1,i+1,iprot))
398 & -atcoord(j,iatsc(1,i,iprot))
400 cacanorm=dsqrt(scalar(caca,caca))
401 c write(iout,*)"i",i," cacanorm",cacanorm
403 caca(k)=caca(k)/cacanorm
408 yy(k)=atcoord(k,iatp(j,i,iprot))-atcoord(k,iatsc(1,i,iprot))
410 xx=scalar(caca,yy)/cacanorm
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))
424 c print *,"i",i," natsc",natsc(i)
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
431 casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot))
435 casc(k)=casc(k)/nnsc-atcoord(k,iatsc(1,i,iprot))
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
445 c write (iout,*) "nnsc",nnsc," iis",iis
446 if (cascnorm.gt.1.0d-2) then
448 casc(k)=casc(k)/cascnorm
454 CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot)
455 & +atforce(k,iatsc(1,i,iprot))
457 do j=2,natsc(i,iprot)
459 if (atsc(j,i,iprot)(1:2).eq."HA") then
462 CGforce_MD(k,i,kkkk,iprot)=CGforce_MD(k,i,kkkk,iprot)
463 & +atforce(k,iatsc(j,i,iprot))
465 c write (4,*) "HA",i,(CGforce_MD(k,i,kkkk),k=1,3)
470 yy(k)=atcoord(k,iatsc(j,i,iprot))-atcoord(k,iatsc(1,i,iprot))
472 xx=scalar(casc,yy)/cascnorm
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))
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)
486 do j=1,natp(nct,iprot)
488 CGforce_MD(k,nct,kkkk,iprot)=CGforce_MD(k,nct,kkkk,iprot)
489 & +atforce(k,iatp(j,nct,iprot))
494 c print *,"nres",nres," nsc",nsc
497 write(iout,*)"CG forces calculated from all-atom forces kkkk",kkkk
498 write (iout,"(4x,a,t55,a)") "Backbone","Sidechain"
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)
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)
514 c-------------------------------------------------------------------
515 subroutine compute_CG_forces_ave(kkkk,iprot)
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
530 c write (iout,'(3i5,3f10.5)') kkkk,j,iatp(j,1,iprot),
531 c & (atforce(k,iatp(j,1,iprot)),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))
540 caca(j)=atcoord(j,iatsc(1,i+1,iprot))
541 & -atcoord(j,iatsc(1,i,iprot))
543 cacanorm=dsqrt(scalar(caca,caca))
544 c write(iout,*)"i",i," cacanorm",cacanorm
546 caca(k)=caca(k)/cacanorm
551 yy(k)=atcoord(k,iatp(j,i,iprot))-atcoord(k,iatsc(1,i,iprot))
553 xx=scalar(caca,yy)/cacanorm
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))
568 c print *,"i",i," natsc",natsc(i)
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
575 casc(k)=casc(k)+atcoord(k,iatsc(j,i,iprot))
579 casc(k)=casc(k)/nnsc-atcoord(k,iatsc(1,i,iprot))
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
589 c write (iout,*) "nnsc",nnsc," iis",iis
590 if (cascnorm.gt.1.0d-2) then
592 casc(k)=casc(k)/cascnorm
598 CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD_ave(k,i,kkkk,iprot)
599 & +atforce(k,iatsc(1,i,iprot))
601 do j=2,natsc(i,iprot)
603 if (atsc(j,i,iprot)(1:2).eq."HA") then
606 CGforce_MD_ave(k,i,kkkk,iprot)=CGforce_MD_ave(k,i,kkkk,iprot)
607 & +atforce(k,iatsc(j,i,iprot))
609 c write (4,*) "HA",i,(CGforce_MD(k,i,kkkk),k=1,3)
614 yy(k)=atcoord(k,iatsc(j,i,iprot))-atcoord(k,iatsc(1,i,iprot))
616 xx=scalar(casc,yy)/cascnorm
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))
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)
631 do j=1,natp(nct,iprot)
633 CGforce_MD_ave(k,nct,kkkk,iprot)=
634 & CGforce_MD_ave(k,nct,kkkk,iprot)
635 & +atforce(k,iatp(j,nct,iprot))
640 c print *,"nres",nres," nsc",nsc
644 c----------------------------------------------------------------------
645 subroutine add_new_MDconf(jjj,jj,jj_old,icount,iprot)
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"
663 include "COMMON.SBRIDGE"
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
669 c write(iout,*) "add_new_MDconf jjj jj jj_old:",jjj,jj,jj_old
671 call int_from_cart1(.false.)
673 c write (iout,*) "After int_from_cart1"
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)
691 & "This conformation WILL NOT be added to the database."
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)
711 & "This conformation WILL NOT be added to the database."
716 if (theta(j).le.0.0d0) then
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)
730 & "This conformation WILL NOT be added to the database."
733 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
737 list_conf_MD(jj,iprot)=jj
738 c write (iout,*) "calling store_MDconf",jjj,jj,icount
740 call store_MDconf_from_file(jj,icount,iprot)
741 c write (iout,*) "finished store_MDconf"
743 if (icount.eq.maxstr_proc) then
745 write (iout,* ) "jj_old",jj_old," jj",jj
746 write (iout,*) "Calling write_and_send_cconf"
749 call write_and_send_MDconf(icount,jj_old,jj,iprot)
752 write (iout,*) "Exited write_and_send_cconf"
759 c------------------------------------------------------------------------------
760 subroutine store_MDconf_from_file(jj,icount,iprot)
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"
771 include "COMMON.FMATCH"
772 integer i,j,jj,icount,ibatch,iprot
773 c Store the conformation that has been read in
776 cg_MD(j,i,icount,iprot)=c(j,i)
779 nss_MD(icount,iprot)=nss
781 ihpb_MD(i,icount,iprot)=ihpb(i)
782 jhpb_MD(i,icount,iprot)=jhpb(i)
784 call compute_CG_forces(icount,iprot)
787 c------------------------------------------------------------------------------
788 subroutine write_and_send_MDconf(icount,jj_old,jj,iprot)
791 include "DIMENSIONS.ZSCOPT"
797 include "COMMON.WEIGHTS"
798 include "COMMON.CHAIN"
799 include "COMMON.SBRIDGE"
800 include "COMMON.INTERACT"
801 include "COMMON.IOUNITS"
802 include "COMMON.CLASSES"
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
814 write (iout,*) "Processor",me," entered WRITE_AND_SEND_MDCONF"
815 write (iout,*) "iprot",iprot
818 call MPI_Bcast(icount,1,MPI_INTEGER,Master,WHAM_COMM,IERROR)
820 write (iout,*) "Processor",me," sent icount=",icount
821 write (iout,*) "Processor",me," jj_old",jj_old," jj",jj
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"
829 call MPI_Bcast(nss_MD(1,iprot),icount,MPI_INTEGER,Master,
831 write (iout,*) "Finished broadcast nss"
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"
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"
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"
848 write (iout,*) "Master finished broadcast"
851 call dawrite_MDcoords(iprot,jj_old,jj,ientout)
853 write (iout,*) "iprot",iprot," jj",jj," nconf_forc",
856 write (iout,*) "Protein",iprot
857 write(iout,*)"Processor",me," received",icount,
858 & " MD conformations and forces"
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),
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,
871 c------------------------------------------------------------------------------
873 subroutine receive_and_write_MDconf(icount,jj_old,jj,iprot)
876 include "DIMENSIONS.ZSCOPT"
878 integer IERROR,STATUS(MPI_STATUS_SIZE)
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"
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
895 do while (icount.gt.0)
897 write (iout,*) "Processor",me," entered RECEIVE_AND_SEND_MDCONF"
898 write (iout,*) "iprot",iprot
901 call MPI_Bcast(icount,1,MPI_INTEGER,Master,WHAM_COMM,IERROR)
903 write (iout,*) "Processor",me," received icount=",icount
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
912 call MPI_Bcast(nss_MD(1,iprot),icount,MPI_INTEGER,Master,
914 write (iout,*) "Finished broadcast nss"
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"
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"
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"
931 write (iout,*) "Processor",me," finished broadcast"
933 call dawrite_MDcoords(iprot,jj_old,jj,ientout)
935 write (iout,*) "Protein",iprot
936 write(iout,*)"Processor",me," received",icount,
937 & " MD conformations and forces"
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),
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,
949 write (iout,*) "jj",jj," nconf_forc",nconf_forc(iprot)
953 c-------------------------------------------------------------------------------
954 subroutine work_partition_MD(lprint)
955 c Split the conformations between processors
958 include "DIMENSIONS.ZSCOPT"
960 include "COMMON.CLASSES"
961 include "COMMON.IOUNITS"
963 include "COMMON.VMCPAR"
964 include "COMMON.CONTROL"
965 integer n,chunk,iprot,i,j,ii,remainder
966 integer kolor,key,ierror,errcode
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.
973 C First try to assign equal number of conformations to each processor.
977 if (.not.fmatch(iprot)) cycle
979 n=ntot_work_MD(iprot)
980 write (iout,*) "Protein",iprot," n=",n
981 indstart_MD(0,iprot)=1
983 scount_MD(0,iprot) = chunk
984 c print *,"i",0," indstart",indstart(0,iprot)," scount",
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",
993 C Determine how many conformations remained yet unassigned.
995 remainder=N-(indstart_MD(nprocs1-1,iprot)
996 & +scount_MD(nprocs1-1,iprot)-1)
997 c print *,"remainder",remainder
999 C Assign the remainder conformations to consecutive processors, starting
1000 C from the lowest rank; this continues until the list is exhausted.
1002 if (remainder .gt. 0) then
1004 scount_MD(i-1,iprot) = scount_MD(i-1,iprot) + 1
1005 indstart_MD(i,iprot) = indstart_MD(i,iprot) + i
1007 do i=remainder+1,nprocs1-1
1008 indstart_MD(i,iprot) = indstart_MD(i,iprot) + remainder
1012 indstart_MD(nprocs1,iprot)=N+1
1013 scount_MD(nprocs1,iprot)=0
1016 indend_MD(i,iprot)=indstart_MD(i,iprot)+scount_MD(i,iprot)-1
1017 idispl_MD(i,iprot)=indstart_MD(i,iprot)-1
1022 N=N+indend_MD(i,iprot)-indstart_MD(i,iprot)+1
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
1029 call MPI_Abort( WHAM_COMM, Ierror, Errcode )
1034 write (iout,*) "Partition of work between processors"
1036 write (iout,*) "Protein",iprot
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)
1047 c------------------------------------------------------------------------------