2 !-----------------------------------------------------------------------
8 !-----------------------------------------------------------------------------
11 !-----------------------------------------------------------------------------
13 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
16 subroutine write_csa_pdb(var,ene,nft,ik,iw_pdb)
19 use geometry_data, only:nres,nvar
20 use geometry, only:var_to_geom,chainbuild
21 use compare, only:secondary2
22 ! implicit real*8 (a-h,o-z)
23 ! include 'DIMENSIONS'
24 ! include 'COMMON.CSA'
25 ! include 'COMMON.BANK'
26 ! include 'COMMON.VAR'
27 ! include 'COMMON.IOUNITS'
28 ! include 'COMMON.MINIM'
29 ! include 'COMMON.SETUP'
30 ! include 'COMMON.GEO'
31 ! include 'COMMON.CHAIN'
32 ! include 'COMMON.LOCAL'
33 ! include 'COMMON.INTERACT'
34 ! include 'COMMON.NAMES'
35 ! include 'COMMON.SBRIDGE'
36 integer :: lenpre,lenpot !,ilen
38 real(kind=8),dimension(nvar) :: var !(maxvar) (maxvar=6*maxres)
39 character(len=50) :: titelloc
40 character(len=3) :: zahl
41 real(kind=8),dimension(mxch*(mxch+1)/2+1) :: ene
43 integer :: nft,ik,iw_pdb
46 if(ene(1).lt.eglob_csa) then
49 call numstr(nglob_csa,zahl)
51 call var_to_geom(nvar,var)
53 call secondary2(.false.)
56 open(icsa_pdb,file=prefix(:lenpre)//'@'//zahl//'.pdb')
59 write(titelloc,'(a2,i3,a3,i9,a3,i6)') &
60 'GM',nglob_csa,' e ',nft,' m ',nmin_csa
62 write(titelloc,'(a2,i3,a3,i9,a3,i6,a5,f5.2,a5,f5.1)') &
63 'GM',nglob_csa,' e ',nft,' m ',nmin_csa,' rms ',&
64 rmsn(ik),' %NC ',pncn(ik)*100
66 call pdbout(eglob_csa,titelloc,icsa_pdb)
71 end subroutine write_csa_pdb
72 !-----------------------------------------------------------------------------
74 !-----------------------------------------------------------------------------
75 subroutine from_pdb(n,idum)
76 ! This subroutine stores the UNRES int variables generated from
77 ! subroutine readpdb into the 1st conformation of in dihang_in.
78 ! Subsequent n-1 conformations of dihang_in have identical values
79 ! of theta and phi as the 1st conformation but random values for
81 ! The array cref (also generated from subroutine readpdb) is stored
82 ! to crefjlee to be used for rmsd calculation in CSA, if necessary.
86 use random, only: ran1
87 ! implicit real*8 (a-h,o-z)
88 ! include 'DIMENSIONS'
89 ! include 'COMMON.IOUNITS'
90 ! include 'COMMON.CHAIN'
91 ! include 'COMMON.VAR'
92 ! include 'COMMON.BANK'
93 ! include 'COMMON.GEO'
95 integer :: n,idum,m,i,j,k,kk,kkk
100 dihang_in(1,j,1,m)=theta(j+1)
101 dihang_in(2,j,1,m)=phi(j+2)
102 dihang_in(3,j,1,m)=alph(j)
103 dihang_in(4,j,1,m)=omeg(j)
105 dihang_in(2,nres-1,1,k)=0.0d0
109 dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
110 dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
111 if(dabs(dihang_in(3,k,1,1)).gt.1.d-6) then
112 dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
113 dihang_in(3,k,1,m)=dihang_in(3,k,1,m)*deg2rad
115 if(dabs(dihang_in(4,k,1,1)).gt.1.d-6) then
116 dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
117 dihang_in(4,k,1,m)=dihang_in(4,k,1,m)*deg2rad
122 ! Store cref to crefjlee (they are in COMMON.CHAIN).
126 crefjlee(kk,k)=cref(kk,k,kkk)
130 open(icsa_native_int,file=csa_native_int,status="old")
132 write(icsa_native_int,*) m,e
133 write(icsa_native_int,200) &
134 (dihang_in(1,k,1,m)*rad2deg,k=2,nres-1)
135 write(icsa_native_int,200) &
136 (dihang_in(2,k,1,m)*rad2deg,k=2,nres-2)
137 write(icsa_native_int,200) &
138 (dihang_in(3,k,1,m)*rad2deg,k=2,nres-1)
139 write(icsa_native_int,200) &
140 (dihang_in(4,k,1,m)*rad2deg,k=2,nres-1)
144 write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
146 close(icsa_native_int)
151 end subroutine from_pdb
152 !-----------------------------------------------------------------------------
153 subroutine from_int(n,mm,idum)
158 use geometry, only:chainbuild,gen_side
159 use energy, only:etotal
161 ! implicit real*8 (a-h,o-z)
162 ! include 'DIMENSIONS'
163 ! include 'COMMON.IOUNITS'
164 ! include 'COMMON.CHAIN'
165 ! include 'COMMON.VAR'
166 ! include 'COMMON.INTERACT'
167 ! include 'COMMON.BANK'
168 ! include 'COMMON.GEO'
169 ! include 'COMMON.CONTACTS'
173 real(kind=8),dimension(0:n_ene) :: energia
175 integer :: n,mm,idum,i,ii,j,m,k,kk,maxcount_fail,icount_fail,maxsi
178 open(icsa_native_int,file=csa_native_int,status="old")
179 read (icsa_native_int,*)
180 call read_angles(icsa_native_int,*10)
182 10 write (iout,'(2a)') "CHUJ NASTAPIL - error in ",&
183 csa_native_int(:ilen(csa_native_int))
187 dihang_in(1,j,1,1)=theta(j+1)
188 dihang_in(2,j,1,1)=phi(j+2)
189 dihang_in(3,j,1,1)=alph(j)
190 dihang_in(4,j,1,1)=omeg(j)
192 dihang_in(2,nres-1,1,1)=0.0d0
194 ! read(icsa_native_int,*) ind,e
195 ! read(icsa_native_int,200) (dihang_in(1,k,1,1),k=2,nres-1)
196 ! read(icsa_native_int,200) (dihang_in(2,k,1,1),k=2,nres-2)
197 ! read(icsa_native_int,200) (dihang_in(3,k,1,1),k=2,nres-1)
198 ! read(icsa_native_int,200) (dihang_in(4,k,1,1),k=2,nres-1)
199 ! dihang_in(2,nres-1,1,1)=0.d0
206 ! dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
207 ! dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
208 ! if(abs(dihang_in(3,k,1,1)).gt.1.d-3) then
209 ! dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
211 ! if(abs(dihang_in(4,k,1,1)).gt.1.d-3) then
212 ! dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
220 DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL)
223 if (itype(i).ne.10) then
224 !d print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
227 do while (fail .and. ii .le. maxsi)
228 call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail)
235 fail = (energia(0).ge.1.0d20)
236 icount_fail=icount_fail+1
240 if (icount_fail.gt.maxcount_fail) then
242 'Failed to generate non-overlaping near-native conf.',&
247 dihang_in(1,j,1,m)=theta(j+1)
248 dihang_in(2,j,1,m)=phi(j+2)
249 dihang_in(3,j,1,m)=alph(j)
250 dihang_in(4,j,1,m)=omeg(j)
252 dihang_in(2,nres-1,1,m)=0.0d0
256 ! write(icsa_native_int,*) m,e
257 ! write(icsa_native_int,200) (dihang_in(1,k,1,m),k=2,nres-1)
258 ! write(icsa_native_int,200) (dihang_in(2,k,1,m),k=2,nres-2)
259 ! write(icsa_native_int,200) (dihang_in(3,k,1,m),k=2,nres-1)
260 ! write(icsa_native_int,200) (dihang_in(4,k,1,m),k=2,nres-1)
262 ! close(icsa_native_int)
267 ! dihang_in(i,j,1,m)=dihang_in(i,j,1,m)*deg2rad
272 call dihang_to_c(dihang_in(1,1,1,1))
274 ! Store c to cref (they are in COMMON.CHAIN).
277 crefjlee(kk,k)=c(kk,k)
281 call contact(.true.,ncont_ref,icont_ref,co)
284 ! write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
286 close(icsa_native_int)
291 end subroutine from_int
292 !-----------------------------------------------------------------------------
293 subroutine dihang_to_c(aarray)
297 use geometry, only:chainbuild
298 ! implicit real*8 (a-h,o-z)
299 ! include 'DIMENSIONS'
300 ! include 'COMMON.CSA'
301 ! include 'COMMON.BANK'
302 ! include 'COMMON.CHAIN'
303 ! include 'COMMON.GEO'
304 ! include 'COMMON.VAR'
306 real(kind=8),dimension(mxang,nres,mxch) :: aarray !(mxang,maxres,mxch)
309 ! phi(i)=dihang_in(1,i-2,1,1)
312 theta(i+1)=aarray(1,i,1)
313 phi(i+2)=aarray(2,i,1)
314 alph(i)=aarray(3,i,1)
315 omeg(i)=aarray(4,i,1)
321 end subroutine dihang_to_c
322 !-----------------------------------------------------------------------------
324 !-----------------------------------------------------------------------------
326 subroutine cartout(time)
328 subroutine cartoutx(time)
330 use geometry_data, only: c,nres
332 use MD_data, only: potE,t_bath
333 ! implicit real*8 (a-h,o-z)
334 ! include 'DIMENSIONS'
335 ! include 'COMMON.CHAIN'
336 ! include 'COMMON.INTERACT'
337 ! include 'COMMON.NAMES'
338 ! include 'COMMON.IOUNITS'
339 ! include 'COMMON.HEADER'
340 ! include 'COMMON.SBRIDGE'
341 ! include 'COMMON.DISTFIT'
342 ! include 'COMMON.MD'
347 #if defined(AIX) || defined(PGI)
348 open(icart,file=cartname,position="append")
350 open(icart,file=cartname,access="append")
352 write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
354 write (icart,'(i4,$)') &
355 nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)
357 write (icart,'(i4,$)') &
358 nss,(ihpb(j),jhpb(j),j=1,nss)
360 write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,&
361 (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),&
362 (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
363 write (icart,'(8f10.5)') &
364 ((c(k,j),k=1,3),j=1,nres),&
365 ((c(k,j+nres),k=1,3),j=nnt,nct)
370 end subroutine cartout
372 end subroutine cartoutx
374 !-----------------------------------------------------------------------------
376 subroutine cartout(time)
377 ! implicit real*8 (a-h,o-z)
378 ! include 'DIMENSIONS'
379 use geometry_data, only: c,nres
381 use MD_data, only: potE,t_bath
385 ! include 'COMMON.SETUP'
387 integer,parameter :: me=0
389 ! include 'COMMON.CHAIN'
390 ! include 'COMMON.INTERACT'
391 ! include 'COMMON.NAMES'
392 ! include 'COMMON.IOUNITS'
393 ! include 'COMMON.HEADER'
394 ! include 'COMMON.SBRIDGE'
395 ! include 'COMMON.DISTFIT'
396 ! include 'COMMON.MD'
400 real(kind=4),dimension(3,2*nres+2) :: xcoord !(3,maxres2+2) (maxres2=2*maxres
405 call xdrfopen_(ixdrf,cartname, "a", iret)
406 call xdrffloat_(ixdrf, real(time), iret)
407 call xdrffloat_(ixdrf, real(potE), iret)
408 call xdrffloat_(ixdrf, real(uconst), iret)
409 call xdrffloat_(ixdrf, real(uconst_back), iret)
410 call xdrffloat_(ixdrf, real(t_bath), iret)
411 call xdrfint_(ixdrf, nss, iret)
414 call xdrfint_(ixdrf, idssb(j)+nres, iret)
415 call xdrfint_(ixdrf, jdssb(j)+nres, iret)
417 call xdrfint_(ixdrf, ihpb(j), iret)
418 call xdrfint_(ixdrf, jhpb(j), iret)
421 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
423 call xdrffloat_(ixdrf, real(qfrag(i)), iret)
426 call xdrffloat_(ixdrf, real(qpair(i)), iret)
429 call xdrffloat_(ixdrf, real(utheta(i)), iret)
430 call xdrffloat_(ixdrf, real(ugamma(i)), iret)
431 call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
434 call xdrfopen(ixdrf,cartname, "a", iret)
435 call xdrffloat(ixdrf, real(time), iret)
436 call xdrffloat(ixdrf, real(potE), iret)
437 call xdrffloat(ixdrf, real(uconst), iret)
438 call xdrffloat(ixdrf, real(uconst_back), iret)
439 call xdrffloat(ixdrf, real(t_bath), iret)
440 call xdrfint(ixdrf, nss, iret)
443 call xdrfint(ixdrf, idssb(j)+nres, iret)
444 call xdrfint(ixdrf, jdssb(j)+nres, iret)
446 call xdrfint(ixdrf, ihpb(j), iret)
447 call xdrfint(ixdrf, jhpb(j), iret)
450 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
452 call xdrffloat(ixdrf, real(qfrag(i)), iret)
455 call xdrffloat(ixdrf, real(qpair(i)), iret)
458 call xdrffloat(ixdrf, real(utheta(i)), iret)
459 call xdrffloat(ixdrf, real(ugamma(i)), iret)
460 call xdrffloat(ixdrf, real(uscdiff(i)), iret)
471 xcoord(j,nres+i-nnt+1)=c(j,i+nres)
477 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
478 call xdrfclose_(ixdrf, iret)
480 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
481 call xdrfclose(ixdrf, iret)
484 end subroutine cartout
486 !-----------------------------------------------------------------------------
487 subroutine statout(itime)
490 use control_data, only:refstr
493 use compare, only:rms_nac_nnc
494 ! implicit real*8 (a-h,o-z)
495 ! include 'DIMENSIONS'
496 ! include 'COMMON.CONTROL'
497 ! include 'COMMON.CHAIN'
498 ! include 'COMMON.INTERACT'
499 ! include 'COMMON.NAMES'
500 ! include 'COMMON.IOUNITS'
501 ! include 'COMMON.HEADER'
502 ! include 'COMMON.SBRIDGE'
503 ! include 'COMMON.DISTFIT'
504 ! include 'COMMON.MD'
505 ! include 'COMMON.REMD'
506 ! include 'COMMON.SETUP'
508 real(kind=8),dimension(0:n_ene) :: energia
509 ! double precision gyrate
511 !el common /gucio/ cm
512 character(len=256) :: line1,line2
513 character(len=4) :: format1,format2
514 character(len=30) :: format
517 real(kind=8) :: rms,frac,frac_nn,co
521 open(istat,file=statname,position="append")
525 open(istat,file=statname,position="append")
527 open(istat,file=statname,access="append")
531 call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
532 write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
533 itime,totT,EK,potE,totE,&
534 rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
537 write (line1,'(i10,f15.2,7f12.3,i5,$)') &
538 itime,totT,EK,potE,totE,&
539 amax,kinetic_T,t_bath,gyrate(),me
542 if(usampl.and.totT.gt.eq_time) then
543 write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,&
544 (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),&
545 (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
546 write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair &
552 if (print_compon) then
554 write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
556 write (istat,format) "#","",&
557 (ename(print_order(i)),i=1,nprint_ene)
559 write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
561 write (istat,format) line1,line2,&
562 (potEcomp(print_order(i)),i=1,nprint_ene)
564 write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
565 write (istat,format) line1,line2
573 end subroutine statout
574 !-----------------------------------------------------------------------------
576 !-----------------------------------------------------------------------------
582 use muca_md, only:read_muca
583 ! implicit real*8 (a-h,o-z)
584 ! include 'DIMENSIONS'
588 ! include 'COMMON.SETUP'
589 ! include 'COMMON.CONTROL'
590 ! include 'COMMON.SBRIDGE'
591 ! include 'COMMON.IOUNITS'
592 logical :: file_exist
594 ! Read force-field parameters except weights
596 ! Read job setup parameters
598 ! Read control parameters for energy minimzation if required
599 if (minim) call read_minim
600 ! Read MCM control parameters if required
601 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
602 ! Read MD control parameters if reqjuired
603 if (modecalc.eq.12) call read_MDpar
604 ! Read MREMD control parameters if required
605 if (modecalc.eq.14) then
609 ! Read MUCA control parameters if required
610 if (lmuca) call read_muca
611 ! Read CSA control parameters if required (from fort.40 if exists
612 ! otherwise from general input file)
613 if (modecalc.eq.8) then
614 inquire (file="fort.40",exist=file_exist)
615 if (.not.file_exist) call csaread
617 !fmc if (modecalc.eq.10) call mcmfread
618 ! Read molecule information, molecule geometry, energy-term weights, and
619 ! restraints if requested
621 ! Print restraint information
623 if (.not. out1file .or. me.eq.king) then
626 write (iout,'(a,i5,a)') "The following",nhpb-nss,&
627 " distance constraints have been imposed"
629 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
634 ! print *,"Processor",myrank," leaves READRTNS"
635 ! write(iout,*) "end readrtns"
637 end subroutine readrtns
638 !-----------------------------------------------------------------------------
641 ! Read molecular data.
643 ! use control, only: ilen
649 use MD_data, only: t_bath
651 use compare, only:seq_comp,contact
653 ! implicit real*8 (a-h,o-z)
654 ! include 'DIMENSIONS'
657 integer :: error_msg,ierror,ierr,ierrcode
659 ! include 'COMMON.IOUNITS'
660 ! include 'COMMON.GEO'
661 ! include 'COMMON.VAR'
662 ! include 'COMMON.INTERACT'
663 ! include 'COMMON.LOCAL'
664 ! include 'COMMON.NAMES'
665 ! include 'COMMON.CHAIN'
666 ! include 'COMMON.FFIELD'
667 ! include 'COMMON.SBRIDGE'
668 ! include 'COMMON.HEADER'
669 ! include 'COMMON.CONTROL'
670 ! include 'COMMON.DBASE'
671 ! include 'COMMON.THREAD'
672 ! include 'COMMON.CONTACTS'
673 ! include 'COMMON.TORCNSTR'
674 ! include 'COMMON.TIME1'
675 ! include 'COMMON.BOUNDS'
676 ! include 'COMMON.MD'
677 ! include 'COMMON.SETUP'
678 character(len=4),dimension(:),allocatable :: sequence !(maxres)
680 ! double precision x(maxvar)
681 character(len=256) :: pdbfile
682 character(len=320) :: weightcard
683 character(len=80) :: weightcard_t!,ucase
684 ! integer,dimension(:),allocatable :: itype_pdb !(maxres)
685 ! common /pizda/ itype_pdb
686 logical :: fail !seq_comp,
687 real(kind=8) :: energia(0:n_ene)
691 integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
693 real(kind=8),dimension(3,maxres2+2) :: c_alloc
694 real(kind=8),dimension(3,0:maxres2) :: dc_alloc
695 integer,dimension(maxres) :: itype_alloc
697 integer :: iti,nsi,maxsi,itrial,itmp
698 real(kind=8) :: wlong,scalscp,co
699 allocate(weights(n_ene))
700 !-----------------------------
701 allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
702 allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
703 allocate(itype(maxres)) !(maxres)
710 !-----------------------------
714 ! Read weights of the subsequent energy terms.
715 call card_concat(weightcard,.true.)
716 call reada(weightcard,'WLONG',wlong,1.0D0)
717 call reada(weightcard,'WSC',wsc,wlong)
718 call reada(weightcard,'WSCP',wscp,wlong)
719 call reada(weightcard,'WELEC',welec,1.0D0)
720 call reada(weightcard,'WVDWPP',wvdwpp,welec)
721 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
722 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
723 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
724 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
725 call reada(weightcard,'WTURN3',wturn3,1.0D0)
726 call reada(weightcard,'WTURN4',wturn4,1.0D0)
727 call reada(weightcard,'WTURN6',wturn6,1.0D0)
728 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
729 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
730 call reada(weightcard,'WBOND',wbond,1.0D0)
731 call reada(weightcard,'WTOR',wtor,1.0D0)
732 call reada(weightcard,'WTORD',wtor_d,1.0D0)
733 call reada(weightcard,'WSHIELD',wshield,0.05D0)
734 call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
735 call reada(weightcard,'WLT',wliptran,1.0D0)
736 call reada(weightcard,'WTUBE',wtube,1.0d0)
737 call reada(weightcard,'WANG',wang,1.0D0)
738 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
739 call reada(weightcard,'SCAL14',scal14,0.4D0)
740 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
741 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
742 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
743 call reada(weightcard,'TEMP0',temp0,300.0d0)
744 if (index(weightcard,'SOFT').gt.0) ipot=6
745 ! 12/1/95 Added weight for the multi-body term WCORR
746 call reada(weightcard,'WCORRH',wcorr,1.0D0)
747 if (wcorr4.gt.0.0d0) wcorr=wcorr4
767 if(me.eq.king.or..not.out1file) &
768 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
769 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
771 10 format (/'Energy-term weights (unscaled):'// &
772 'WSCC= ',f10.6,' (SC-SC)'/ &
773 'WSCP= ',f10.6,' (SC-p)'/ &
774 'WELEC= ',f10.6,' (p-p electr)'/ &
775 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
776 'WBOND= ',f10.6,' (stretching)'/ &
777 'WANG= ',f10.6,' (bending)'/ &
778 'WSCLOC= ',f10.6,' (SC local)'/ &
779 'WTOR= ',f10.6,' (torsional)'/ &
780 'WTORD= ',f10.6,' (double torsional)'/ &
781 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
782 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
783 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
784 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
785 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
786 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
787 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
788 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
789 'WTURN6= ',f10.6,' (turns, 6th order)')
790 if(me.eq.king.or..not.out1file)then
791 if (wcorr4.gt.0.0d0) then
792 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
793 'between contact pairs of peptide groups'
794 write (iout,'(2(a,f5.3/))') &
795 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
796 'Range of quenching the correlation terms:',2*delt_corr
797 else if (wcorr.gt.0.0d0) then
798 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
799 'between contact pairs of peptide groups'
801 write (iout,'(a,f8.3)') &
802 'Scaling factor of 1,4 SC-p interactions:',scal14
803 write (iout,'(a,f8.3)') &
804 'General scaling factor of SC-p interactions:',scalscp
806 r0_corr=cutoff_corr-delt_corr
808 aad(i,1)=scalscp*aad(i,1)
809 aad(i,2)=scalscp*aad(i,2)
810 bad(i,1)=scalscp*bad(i,1)
811 bad(i,2)=scalscp*bad(i,2)
813 call rescale_weights(t_bath)
814 if(me.eq.king.or..not.out1file) &
815 write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
816 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
818 22 format (/'Energy-term weights (scaled):'// &
819 'WSCC= ',f10.6,' (SC-SC)'/ &
820 'WSCP= ',f10.6,' (SC-p)'/ &
821 'WELEC= ',f10.6,' (p-p electr)'/ &
822 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
823 'WBOND= ',f10.6,' (stretching)'/ &
824 'WANG= ',f10.6,' (bending)'/ &
825 'WSCLOC= ',f10.6,' (SC local)'/ &
826 'WTOR= ',f10.6,' (torsional)'/ &
827 'WTORD= ',f10.6,' (double torsional)'/ &
828 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
829 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
830 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
831 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
832 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
833 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
834 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
835 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
836 'WTURN6= ',f10.6,' (turns, 6th order)')
837 if(me.eq.king.or..not.out1file) &
838 write (iout,*) "Reference temperature for weights calculation:",&
840 call reada(weightcard,"D0CM",d0cm,3.78d0)
841 call reada(weightcard,"AKCM",akcm,15.1d0)
842 call reada(weightcard,"AKTH",akth,11.0d0)
843 call reada(weightcard,"AKCT",akct,12.0d0)
844 call reada(weightcard,"V1SS",v1ss,-1.08d0)
845 call reada(weightcard,"V2SS",v2ss,7.61d0)
846 call reada(weightcard,"V3SS",v3ss,13.7d0)
847 call reada(weightcard,"EBR",ebr,-5.50D0)
848 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
850 call reada(weightcard,"HT",Ht,0.0D0)
852 ss_depth=ebr/wsc-0.25*eps(1,1)
853 Ht=Ht/wsc-0.25*eps(1,1)
854 akcm=akcm*wstrain/wsc
855 akth=akth*wstrain/wsc
856 akct=akct*wstrain/wsc
857 v1ss=v1ss*wstrain/wsc
858 v2ss=v2ss*wstrain/wsc
859 v3ss=v3ss*wstrain/wsc
861 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
864 if(me.eq.king.or..not.out1file) then
865 write (iout,*) "Parameters of the SS-bond potential:"
866 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
868 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
869 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
870 write (iout,*)" HT",Ht
871 print *,'indpdb=',indpdb,' pdbref=',pdbref
873 if (indpdb.gt.0 .or. pdbref) then
874 read(inp,'(a)') pdbfile
875 if(me.eq.king.or..not.out1file) &
876 write (iout,'(2a)') 'PDB data will be read from file ',&
877 pdbfile(:ilen(pdbfile))
878 open(ipdbin,file=pdbfile,status='old',err=33)
880 33 write (iout,'(a)') 'Error opening PDB file.'
883 ! print *,'Begin reading pdb data'
885 ! print *,'Finished reading pdb data'
886 if(me.eq.king.or..not.out1file) &
887 write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
888 ' nstart_sup=',nstart_sup !,"ergwergewrgae"
889 !el if(.not.allocated(itype_pdb))
890 allocate(itype_pdb(nres))
892 itype_pdb(i)=itype(i)
896 nct=nstart_sup+nsup-1
897 !el if(.not.allocated(icont_ref))
898 allocate(icont_ref(2,12*nres)) ! maxcont=12*maxres
899 call contact(.false.,ncont_ref,icont_ref,co)
902 if(me.eq.king.or..not.out1file) &
903 write(iout,*)'Adding sidechains'
907 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
910 do while (fail.and.nsi.le.maxsi)
911 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
914 if(fail) write(iout,*)'Adding sidechain failed for res ',&
915 i,' after ',nsi,' trials'
921 if (indpdb.eq.0) then
922 ! Read sequence if not taken from the pdb file.
924 ! print *,'nres=',nres
925 allocate(sequence(nres))
926 if (iscode.gt.0) then
927 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
929 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
931 ! Convert sequence to numeric code
933 itype(i)=rescode(i,sequence(i),iscode)
935 ! Assign initial virtual bond lengths
936 !elwrite(iout,*) "test_alloc"
937 if(.not.allocated(vbld)) allocate(vbld(2*nres))
938 !elwrite(iout,*) "test_alloc"
939 if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
940 !elwrite(iout,*) "test_alloc"
946 vbld(i+nres)=dsc(iabs(itype(i)))
947 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
948 ! write (iout,*) "i",i," itype",itype(i),
949 ! & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
953 ! print '(20i4)',(itype(i),i=1,nres)
954 !----------------------------
955 !el reallocate tables
958 ! c_alloc(j,i)=c(j,i)
959 ! dc_alloc(j,i)=dc(j,i)
963 !elwrite(iout,*) "itype",i,itype(i)
964 ! itype_alloc(i)=itype(i)
970 ! allocate(c(3,2*nres+4))
971 ! allocate(dc(3,0:2*nres+2))
972 ! allocate(itype(nres+2))
973 allocate(itel(nres+2))
978 ! c(j,i)=c_alloc(j,i)
979 ! dc(j,i)=dc_alloc(j,i)
983 ! itype(i)=itype_alloc(i)
986 !--------------------------
989 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
991 if (itype(i).eq.ntyp1) then
995 else if (iabs(itype(i+1)).ne.20) then
997 else if (iabs(itype(i)).ne.20) then
1004 if(me.eq.king.or..not.out1file)then
1005 write (iout,*) "ITEL"
1007 write (iout,*) i,itype(i),itel(i)
1009 print *,'Call Read_Bridge.'
1012 !--------------------------------
1013 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1014 call alloc_geo_arrays
1015 call alloc_ener_arrays
1016 !--------------------------------
1017 ! 8/13/98 Set limits to generating the dihedral angles
1022 read (inp,*) ndih_constr
1023 if (ndih_constr.gt.0) then
1024 allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1025 allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1027 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1028 if(me.eq.king.or..not.out1file)then
1030 'There are',ndih_constr,' constraints on phi angles.'
1032 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1036 phi0(i)=deg2rad*phi0(i)
1037 drange(i)=deg2rad*drange(i)
1039 if(me.eq.king.or..not.out1file) &
1040 write (iout,*) 'FTORS',ftors
1043 phibound(1,ii) = phi0(i)-drange(i)
1044 phibound(2,ii) = phi0(i)+drange(i)
1049 if (me.eq.king) then
1051 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1053 write (iout,'(a3,i5,2f10.1)') &
1054 restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1060 !d print *,'NNT=',NNT,' NCT=',NCT
1061 if (itype(1).eq.ntyp1) nnt=2
1062 if (itype(nres).eq.ntyp1) nct=nct-1
1064 if(me.eq.king.or..not.out1file) &
1065 write (iout,'(a,i3)') 'nsup=',nsup
1067 if (nsup.le.(nct-nnt+1)) then
1068 do i=0,nct-nnt+1-nsup
1069 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1074 write (iout,'(a)') &
1075 'Error - sequences to be superposed do not match.'
1078 do i=0,nsup-(nct-nnt+1)
1079 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1081 nstart_sup=nstart_sup+i
1086 write (iout,'(a)') &
1087 'Error - sequences to be superposed do not match.'
1090 if (nsup.eq.0) nsup=nct-nnt
1091 if (nstart_sup.eq.0) nstart_sup=nnt
1092 if (nstart_seq.eq.0) nstart_seq=nnt
1093 if(me.eq.king.or..not.out1file) &
1094 write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1095 ' nstart_seq=',nstart_seq !,"242343453254"
1097 !--- Zscore rms -------
1098 if (nz_start.eq.0) nz_start=nnt
1099 if (nz_end.eq.0 .and. nsup.gt.0) then
1101 else if (nz_end.eq.0) then
1104 if(me.eq.king.or..not.out1file)then
1105 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1106 write (iout,*) 'IZ_SC=',iz_sc
1108 !----------------------
1111 if (.not.pdbref) then
1112 call read_angles(inp,*38)
1114 38 write (iout,'(a)') 'Error reading reference structure.'
1116 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1117 stop 'Error reading reference structure'
1121 !zscore call geom_to_var(nvar,coord_exp_zs(1,1))
1128 cref(j,i,kkk)=c(j,i)
1131 call contact(.true.,ncont_ref,icont_ref,co)
1133 ! write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1135 if (constr_dist.gt.0) call read_dist_constr
1136 write (iout,*) "After read_dist_constr nhpb",nhpb
1138 if(me.eq.king.or..not.out1file) &
1139 write (iout,*) 'Contact order:',co
1141 if(me.eq.king.or..not.out1file) &
1142 write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1145 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1147 if(me.eq.king.or..not.out1file) &
1148 write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',&
1149 icont_ref(1,i),' ',&
1150 restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1154 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1155 .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1156 modecalc.ne.10) then
1157 ! If input structure hasn't been supplied from the PDB file read or generate
1159 if (iranconf.eq.0 .and. .not. extconf) then
1160 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1161 write (iout,'(a)') 'Initial geometry will be read in.'
1163 read(inp,'(8f10.5)',end=36,err=36) &
1164 ((c(l,k),l=1,3),k=1,nres),&
1165 ((c(l,k+nres),l=1,3),k=nnt,nct)
1166 write (iout,*) "Exit READ_CART"
1167 write (iout,'(8f10.5)') &
1168 ((c(l,k),l=1,3),k=1,nres),&
1169 ((c(l,k+nres),l=1,3),k=nnt,nct)
1170 call int_from_cart1(.true.)
1171 write (iout,*) "Finish INT_TO_CART"
1174 dc(j,i)=c(j,i+1)-c(j,i)
1175 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1179 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1181 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1182 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1188 call read_angles(inp,*36)
1191 36 write (iout,'(a)') 'Error reading angle file.'
1193 call mpi_finalize( MPI_COMM_WORLD,IERR )
1195 stop 'Error reading angle file.'
1197 else if (extconf) then
1198 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1199 write (iout,'(a)') 'Extended chain initial geometry.'
1201 theta(i)=90d0*deg2rad
1204 phi(i)=180d0*deg2rad
1207 alph(i)=110d0*deg2rad
1209 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1211 omeg(i)=-120d0*deg2rad
1212 if (itype(i).le.0) omeg(i)=-omeg(i)
1215 if(me.eq.king.or..not.out1file) &
1216 write (iout,'(a)') 'Random-generated initial geometry.'
1220 if (me.eq.king .or. fg_rank.eq.0 .and. &
1221 ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1225 call gen_rand_conf(itmp,*30)
1227 30 write (iout,*) 'Failed to generate random conformation',&
1229 write (*,*) 'Processor:',me,&
1230 ' Failed to generate random conformation',&
1240 write (iout,'(a,i3,a)') 'Processor:',me,&
1241 ' error in generating random conformation.'
1242 write (*,'(a,i3,a)') 'Processor:',me,&
1243 ' error in generating random conformation.'
1246 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1252 call gen_rand_conf(itmp,*335)
1254 335 write (iout,*) 'Failed to generate random conformation',&
1256 write (*,*) 'Failed to generate random conformation',&
1259 write (iout,'(a,i3,a)') 'Processor:',me,&
1260 ' error in generating random conformation.'
1261 write (*,'(a,i3,a)') 'Processor:',me,&
1262 ' error in generating random conformation.'
1267 elseif (modecalc.eq.4) then
1268 read (inp,'(a)') intinname
1269 open (intin,file=intinname,status='old',err=333)
1270 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1271 write (iout,'(a)') 'intinname',intinname
1272 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1274 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1276 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1278 stop 'Error opening angle file.'
1282 ! Generate distance constraints, if the PDB structure is to be regularized.
1283 if (nthread.gt.0) then
1284 call read_threadbase
1287 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1288 if (me.eq.king .or. .not. out1file) &
1290 !elwrite (iout,*)"alph(i)*rad2deg",(alph(i)*rad2deg, i=1,nres)
1291 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1292 write (iout,'(/a,i3,a)') &
1293 'The chain contains',ns,' disulfide-bridging cysteines.'
1294 write (iout,'(20i4)') (iss(i),i=1,ns)
1296 write(iout,*)"Running with dynamic disulfide-bond formation"
1298 write (iout,'(/a/)') 'Pre-formed links are:'
1304 if (me.eq.king.or..not.out1file) &
1305 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1306 restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),&
1312 if (ns.gt.0.and.dyn_ss) then
1316 forcon(i-nss)=forcon(i)
1323 dyn_ss_mask(iss(i))=.true.
1326 if (i2ndstr.gt.0) call secstrp2dihc
1327 ! call geom_to_var(nvar,x)
1328 ! call etotal(energia(0))
1329 ! call enerprint(energia(0))
1330 ! call briefout(0,etot)
1332 !d write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1333 !d write (iout,'(a)') 'Variable list:'
1334 !d write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1336 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1337 write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1338 'Processor',myrank,': end reading molecular data.'
1341 end subroutine molread
1342 !-----------------------------------------------------------------------------
1343 !-----------------------------------------------------------------------------