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,1).ne.10) then
224 !d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
227 do while (fail .and. ii .le. maxsi)
228 call gen_side(itype(i,1),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)
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
516 integer :: i,ii1,ii2,j
517 real(kind=8) :: rms,frac,frac_nn,co,distance
521 open(istat,file=statname,position="append")
525 open(istat,file=statname,position="append")
527 open(istat,file=statname,access="append")
530 if (AFMlog.gt.0) then
532 call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
533 write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')&
534 itime,totT,EK,potE,totE,&
535 rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
539 !C print *,'A CHUJ',potEcomp(23)
540 write (line1,'(i10,f15.2,7f12.3,i5,$)') &
541 itime,totT,EK,potE,totE,&
542 kinetic_T,t_bath,gyrate(),&
546 else if (selfguide.gt.0) then
549 distance=distance+(c(j,afmend)-c(j,afmbeg))**2
551 distance=dsqrt(distance)
553 call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
554 write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
556 itime,totT,EK,potE,totE,&
557 rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
558 distance,potEcomp(23),me
562 !C print *,'A CHUJ',potEcomp(23)
563 write (line1,'(i10,f15.2,8f12.3,i5,$)')&
564 itime,totT,EK,potE,totE, &
565 kinetic_T,t_bath,gyrate(),&
566 distance,potEcomp(23),me
571 call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
572 write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
573 itime,totT,EK,potE,totE,&
574 rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
577 write (line1,'(i10,f15.2,7f12.3,i5,$)') &
578 itime,totT,EK,potE,totE,&
579 amax,kinetic_T,t_bath,gyrate(),me
583 if(usampl.and.totT.gt.eq_time) then
584 write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,&
585 (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),&
586 (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
587 write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair &
593 if (print_compon) then
595 write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
597 write (istat,format) "#","",&
598 (ename(print_order(i)),i=1,nprint_ene)
600 write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
602 write (istat,format) line1,line2,&
603 (potEcomp(print_order(i)),i=1,nprint_ene)
605 write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
606 write (istat,format) line1,line2
614 end subroutine statout
615 !-----------------------------------------------------------------------------
617 !-----------------------------------------------------------------------------
623 use muca_md, only:read_muca
624 ! implicit real*8 (a-h,o-z)
625 ! include 'DIMENSIONS'
629 ! include 'COMMON.SETUP'
630 ! include 'COMMON.CONTROL'
631 ! include 'COMMON.SBRIDGE'
632 ! include 'COMMON.IOUNITS'
633 logical :: file_exist
635 ! Read force-field parameters except weights
637 ! Read job setup parameters
639 ! Read control parameters for energy minimzation if required
640 if (minim) call read_minim
641 ! Read MCM control parameters if required
642 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
643 ! Read MD control parameters if reqjuired
644 if (modecalc.eq.12) call read_MDpar
645 ! Read MREMD control parameters if required
646 if (modecalc.eq.14) then
650 ! Read MUCA control parameters if required
651 if (lmuca) call read_muca
652 ! Read CSA control parameters if required (from fort.40 if exists
653 ! otherwise from general input file)
654 if (modecalc.eq.8) then
655 inquire (file="fort.40",exist=file_exist)
656 if (.not.file_exist) call csaread
658 !fmc if (modecalc.eq.10) call mcmfread
659 ! Read molecule information, molecule geometry, energy-term weights, and
660 ! restraints if requested
662 ! Print restraint information
664 if (.not. out1file .or. me.eq.king) then
667 write (iout,'(a,i5,a)') "The following",nhpb-nss,&
668 " distance constraints have been imposed"
670 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
675 ! print *,"Processor",myrank," leaves READRTNS"
676 ! write(iout,*) "end readrtns"
678 end subroutine readrtns
679 !-----------------------------------------------------------------------------
682 ! Read molecular data.
684 ! use control, only: ilen
690 use MD_data, only: t_bath
692 use compare, only:seq_comp,contact
694 ! implicit real*8 (a-h,o-z)
695 ! include 'DIMENSIONS'
698 integer :: error_msg,ierror,ierr,ierrcode
700 ! include 'COMMON.IOUNITS'
701 ! include 'COMMON.GEO'
702 ! include 'COMMON.VAR'
703 ! include 'COMMON.INTERACT'
704 ! include 'COMMON.LOCAL'
705 ! include 'COMMON.NAMES'
706 ! include 'COMMON.CHAIN'
707 ! include 'COMMON.FFIELD'
708 ! include 'COMMON.SBRIDGE'
709 ! include 'COMMON.HEADER'
710 ! include 'COMMON.CONTROL'
711 ! include 'COMMON.DBASE'
712 ! include 'COMMON.THREAD'
713 ! include 'COMMON.CONTACTS'
714 ! include 'COMMON.TORCNSTR'
715 ! include 'COMMON.TIME1'
716 ! include 'COMMON.BOUNDS'
717 ! include 'COMMON.MD'
718 ! include 'COMMON.SETUP'
719 character(len=4),dimension(:,:),allocatable :: sequence !(maxres,maxmolecules)
721 ! double precision x(maxvar)
722 character(len=256) :: pdbfile
723 character(len=800) :: weightcard
724 character(len=80) :: weightcard_t!,ucase
725 ! integer,dimension(:),allocatable :: itype_pdb !(maxres)
726 ! common /pizda/ itype_pdb
727 logical :: fail !seq_comp,
728 real(kind=8) :: energia(0:n_ene)
732 integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
734 real(kind=8),dimension(3,maxres2+2) :: c_alloc
735 real(kind=8),dimension(3,0:maxres2) :: dc_alloc
736 integer,dimension(maxres) :: itype_alloc
738 integer :: iti,nsi,maxsi,itrial,itmp
739 real(kind=8) :: wlong,scalscp,co,ssscale
740 allocate(weights(n_ene))
741 !-----------------------------
742 allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
743 allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
744 allocate(itype(maxres,5)) !(maxres)
745 allocate(istype(maxres))
752 !-----------------------------
756 ! Read weights of the subsequent energy terms.
757 call card_concat(weightcard,.true.)
758 call reada(weightcard,'WLONG',wlong,1.0D0)
759 call reada(weightcard,'WSC',wsc,wlong)
760 call reada(weightcard,'WSCP',wscp,wlong)
761 call reada(weightcard,'WELEC',welec,1.0D0)
762 call reada(weightcard,'WVDWPP',wvdwpp,welec)
763 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
764 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
765 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
766 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
767 call reada(weightcard,'WTURN3',wturn3,1.0D0)
768 call reada(weightcard,'WTURN4',wturn4,1.0D0)
769 call reada(weightcard,'WTURN6',wturn6,1.0D0)
770 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
771 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
772 call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,1.0D0)
773 call reada(weightcard,'WELPP',welpp,1.0d0)
774 call reada(weightcard,'WVDWPSB',wvdwpsb,1.0d0)
775 call reada(weightcard,'WELPSB',welpsb,1.0D0)
776 call reada(weightcard,'WVDWSB',wvdwsb,1.0d0)
777 call reada(weightcard,'WELSB',welsb,1.0D0)
778 call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0)
779 call reada(weightcard,'WANG_NUCL',wang_nucl,1.0D0)
780 call reada(weightcard,'WSBLOC',wsbloc,1.0D0)
781 call reada(weightcard,'WTOR_NUCL',wtor_nucl,1.0D0)
782 call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,1.0D0)
783 call reada(weightcard,'WCORR_NUCL',wcorr_nucl,1.0D0)
784 call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,1.0D0)
785 call reada(weightcard,'WBOND',wbond,1.0D0)
786 call reada(weightcard,'WTOR',wtor,1.0D0)
787 call reada(weightcard,'WTORD',wtor_d,1.0D0)
788 call reada(weightcard,'WSHIELD',wshield,0.05D0)
789 call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
790 call reada(weightcard,'WLT',wliptran,1.0D0)
791 call reada(weightcard,'WTUBE',wtube,1.0d0)
792 call reada(weightcard,'WANG',wang,1.0D0)
793 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
794 call reada(weightcard,'SCAL14',scal14,0.4D0)
795 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
796 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
797 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
798 call reada(weightcard,'TEMP0',temp0,300.0d0)
799 call reada(weightcard,'WSCBASE',wscbase,1.0D0)
800 if (index(weightcard,'SOFT').gt.0) ipot=6
801 call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0)
802 call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
803 call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
804 call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
805 call reada(weightcard,'WSCPHO',wscpho,1.0d0)
806 call reada(weightcard,'WPEPPHO',wpeppho,1.0d0)
808 ! 12/1/95 Added weight for the multi-body term WCORR
809 call reada(weightcard,'WCORRH',wcorr,1.0D0)
810 if (wcorr4.gt.0.0d0) wcorr=wcorr4
830 if(me.eq.king.or..not.out1file) &
831 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
832 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
834 10 format (/'Energy-term weights (unscaled):'// &
835 'WSCC= ',f10.6,' (SC-SC)'/ &
836 'WSCP= ',f10.6,' (SC-p)'/ &
837 'WELEC= ',f10.6,' (p-p electr)'/ &
838 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
839 'WBOND= ',f10.6,' (stretching)'/ &
840 'WANG= ',f10.6,' (bending)'/ &
841 'WSCLOC= ',f10.6,' (SC local)'/ &
842 'WTOR= ',f10.6,' (torsional)'/ &
843 'WTORD= ',f10.6,' (double torsional)'/ &
844 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
845 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
846 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
847 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
848 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
849 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
850 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
851 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
852 'WTURN6= ',f10.6,' (turns, 6th order)')
853 if(me.eq.king.or..not.out1file)then
854 if (wcorr4.gt.0.0d0) then
855 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
856 'between contact pairs of peptide groups'
857 write (iout,'(2(a,f5.3/))') &
858 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
859 'Range of quenching the correlation terms:',2*delt_corr
860 else if (wcorr.gt.0.0d0) then
861 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
862 'between contact pairs of peptide groups'
864 write (iout,'(a,f8.3)') &
865 'Scaling factor of 1,4 SC-p interactions:',scal14
866 write (iout,'(a,f8.3)') &
867 'General scaling factor of SC-p interactions:',scalscp
869 r0_corr=cutoff_corr-delt_corr
871 aad(i,1)=scalscp*aad(i,1)
872 aad(i,2)=scalscp*aad(i,2)
873 bad(i,1)=scalscp*bad(i,1)
874 bad(i,2)=scalscp*bad(i,2)
876 call rescale_weights(t_bath)
877 if(me.eq.king.or..not.out1file) &
878 write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
879 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
881 22 format (/'Energy-term weights (scaled):'// &
882 'WSCC= ',f10.6,' (SC-SC)'/ &
883 'WSCP= ',f10.6,' (SC-p)'/ &
884 'WELEC= ',f10.6,' (p-p electr)'/ &
885 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
886 'WBOND= ',f10.6,' (stretching)'/ &
887 'WANG= ',f10.6,' (bending)'/ &
888 'WSCLOC= ',f10.6,' (SC local)'/ &
889 'WTOR= ',f10.6,' (torsional)'/ &
890 'WTORD= ',f10.6,' (double torsional)'/ &
891 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
892 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
893 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
894 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
895 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
896 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
897 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
898 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
899 'WTURN6= ',f10.6,' (turns, 6th order)')
900 if(me.eq.king.or..not.out1file) &
901 write (iout,*) "Reference temperature for weights calculation:",&
903 call reada(weightcard,"D0CM",d0cm,3.78d0)
904 call reada(weightcard,"AKCM",akcm,15.1d0)
905 call reada(weightcard,"AKTH",akth,11.0d0)
906 call reada(weightcard,"AKCT",akct,12.0d0)
907 call reada(weightcard,"V1SS",v1ss,-1.08d0)
908 call reada(weightcard,"V2SS",v2ss,7.61d0)
909 call reada(weightcard,"V3SS",v3ss,13.7d0)
910 call reada(weightcard,"EBR",ebr,-5.50D0)
911 call reada(weightcard,"ATRISS",atriss,0.301D0)
912 call reada(weightcard,"BTRISS",btriss,0.021D0)
913 call reada(weightcard,"CTRISS",ctriss,1.001D0)
914 call reada(weightcard,"DTRISS",dtriss,1.001D0)
915 call reada(weightcard,"SSSCALE",ssscale,1.0D0)
916 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
918 call reada(weightcard,"HT",Ht,0.0D0)
920 ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
921 Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
922 akcm=akcm/wsc*ssscale
923 akth=akth/wsc*ssscale
924 akct=akct/wsc*ssscale
925 v1ss=v1ss/wsc*ssscale
926 v2ss=v2ss/wsc*ssscale
927 v3ss=v3ss/wsc*ssscale
929 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
932 if(me.eq.king.or..not.out1file) then
933 write (iout,*) "Parameters of the SS-bond potential:"
934 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
936 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
937 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
938 write (iout,*)" HT",Ht
939 print *,'indpdb=',indpdb,' pdbref=',pdbref
941 if (indpdb.gt.0 .or. pdbref) then
942 read(inp,'(a)') pdbfile
943 if(me.eq.king.or..not.out1file) &
944 write (iout,'(2a)') 'PDB data will be read from file ',&
945 pdbfile(:ilen(pdbfile))
946 open(ipdbin,file=pdbfile,status='old',err=33)
948 33 write (iout,'(a)') 'Error opening PDB file.'
951 ! print *,'Begin reading pdb data'
953 ! print *,'Finished reading pdb data'
954 if(me.eq.king.or..not.out1file) &
955 write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
956 ' nstart_sup=',nstart_sup !,"ergwergewrgae"
957 !el if(.not.allocated(itype_pdb))
958 allocate(itype_pdb(nres))
960 itype_pdb(i)=itype(i,1)
964 nct=nstart_sup+nsup-1
965 !el if(.not.allocated(icont_ref))
966 allocate(icont_ref(2,(nres/2)*nres)) ! maxcont=12*maxres
967 call contact(.false.,ncont_ref,icont_ref,co)
970 if(me.eq.king.or..not.out1file) &
971 write(iout,*)'Adding sidechains'
975 if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
978 do while (fail.and.nsi.le.maxsi)
979 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
982 if(fail) write(iout,*)'Adding sidechain failed for res ',&
983 i,' after ',nsi,' trials'
989 if (indpdb.eq.0) then
991 allocate(sequence(maxres,5))
995 ! Read sequence if not taken from the pdb file.
997 read (inp,*) nres_molec(molec)
999 if (iscode.gt.0) then
1000 read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
1002 read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1004 ! read(inp,*) weightcard_t
1005 ! print *,"po seq" weightcard_t
1006 ! Convert sequence to numeric code
1008 do i=1,nres_molec(molec)
1010 itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1015 ! read(inp,*) weightcard_t
1016 ! print *,"po seq", weightcard_t
1019 ! Read sequence if not taken from the pdb file.
1021 read (inp,*) nres_molec(molec)
1022 ! print *,'nres=',nres
1023 ! allocate(sequence(maxres,5))
1024 ! if (iscode.gt.0) then
1025 read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
1026 ! Convert sequence to numeric code
1028 do i=1,nres_molec(molec)
1030 istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)
1031 itype(itmp,molec)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
1036 ! Read sequence if not taken from the pdb file.
1038 read (inp,*) nres_molec(molec)
1039 ! print *,'nres=',nres
1040 read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1041 ! Convert sequence to numeric code
1042 print *,nres_molec(molec)
1043 do i=1,nres_molec(molec)
1046 itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1051 nres=nres+nres_molec(i)
1052 print *,nres_molec(i)
1055 ! Assign initial virtual bond lengths
1056 if(.not.allocated(molnum)) then
1057 allocate(molnum(nres+1))
1060 do j=1,nres_molec(i)
1065 ! print *,nres_molec(i)
1067 if(.not.allocated(vbld)) allocate(vbld(2*nres))
1068 if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1074 print *, "molnum",molnum(i),itype(i,molnum(i))
1075 vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
1076 vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
1077 ! write (iout,*) "i",i," itype",itype(i,1),
1078 ! & " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1082 ! print '(20i4)',(itype(i,1),i=1,nres)
1083 !----------------------------
1084 !el reallocate tables
1087 ! c_alloc(j,i)=c(j,i)
1088 ! dc_alloc(j,i)=dc(j,i)
1092 !elwrite(iout,*) "itype",i,itype(i,1)
1093 ! itype_alloc(i)=itype(i,1)
1099 ! allocate(c(3,2*nres+4))
1100 ! allocate(dc(3,0:2*nres+2))
1101 ! allocate(itype(nres+2))
1102 allocate(itel(nres+2))
1107 ! c(j,i)=c_alloc(j,i)
1108 ! dc(j,i)=dc_alloc(j,i)
1112 ! itype(i,1)=itype_alloc(i)
1115 !--------------------------
1118 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1120 if (itype(i,1).eq.ntyp1) then
1124 else if (iabs(itype(i+1,1)).ne.20) then
1126 else if (iabs(itype(i,1)).ne.20) then
1133 if(me.eq.king.or..not.out1file)then
1134 write (iout,*) "ITEL"
1137 write (iout,*) i,itype(i,1),itel(i)
1139 print *,'Call Read_Bridge.'
1142 !--------------------------------
1143 print *,"tu dochodze"
1144 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1145 call alloc_geo_arrays
1146 call alloc_ener_arrays
1147 !--------------------------------
1148 ! 8/13/98 Set limits to generating the dihedral angles
1153 read (inp,*) ndih_constr
1154 if (ndih_constr.gt.0) then
1155 allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1156 allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1159 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1160 if(me.eq.king.or..not.out1file)then
1162 'There are',ndih_constr,' constraints on phi angles.'
1164 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1168 phi0(i)=deg2rad*phi0(i)
1169 drange(i)=deg2rad*drange(i)
1171 if(me.eq.king.or..not.out1file) &
1172 write (iout,*) 'FTORS',ftors
1175 phibound(1,ii) = phi0(i)-drange(i)
1176 phibound(2,ii) = phi0(i)+drange(i)
1179 if (with_theta_constr) then
1180 !C with_theta_constr is keyword allowing for occurance of theta constrains
1181 read (inp,*) ntheta_constr
1182 !C ntheta_constr is the number of theta constrains
1183 if (ntheta_constr.gt.0) then
1184 !C read (inp,*) ftors
1185 allocate(itheta_constr(ntheta_constr))
1186 allocate(theta_constr0(ntheta_constr))
1187 allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1188 read (inp,*) (itheta_constr(i),theta_constr0(i), &
1189 theta_drange(i),for_thet_constr(i), &
1191 !C the above code reads from 1 to ntheta_constr
1192 !C itheta_constr(i) residue i for which is theta_constr
1193 !C theta_constr0 the global minimum value
1194 !C theta_drange is range for which there is no energy penalty
1195 !C for_thet_constr is the force constant for quartic energy penalty
1197 if(me.eq.king.or..not.out1file)then
1199 'There are',ntheta_constr,' constraints on phi angles.'
1200 do i=1,ntheta_constr
1201 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1206 do i=1,ntheta_constr
1207 theta_constr0(i)=deg2rad*theta_constr0(i)
1208 theta_drange(i)=deg2rad*theta_drange(i)
1210 !C if(me.eq.king.or..not.out1file)
1211 !C & write (iout,*) 'FTORS',ftors
1212 !C do i=1,ntheta_constr
1213 !C ii = itheta_constr(i)
1214 !C thetabound(1,ii) = phi0(i)-drange(i)
1215 !C thetabound(2,ii) = phi0(i)+drange(i)
1217 endif ! ntheta_constr.gt.0
1218 endif! with_theta_constr
1222 if (me.eq.king) then
1224 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1226 write (iout,'(a3,i5,2f10.1)') &
1227 restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1233 print *,'NNT=',NNT,' NCT=',NCT
1234 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1235 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1237 if(me.eq.king.or..not.out1file) &
1238 write (iout,'(a,i3)') 'nsup=',nsup
1240 if (nsup.le.(nct-nnt+1)) then
1241 do i=0,nct-nnt+1-nsup
1242 if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1247 write (iout,'(a)') &
1248 'Error - sequences to be superposed do not match.'
1251 do i=0,nsup-(nct-nnt+1)
1252 if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1254 nstart_sup=nstart_sup+i
1259 write (iout,'(a)') &
1260 'Error - sequences to be superposed do not match.'
1263 if (nsup.eq.0) nsup=nct-nnt
1264 if (nstart_sup.eq.0) nstart_sup=nnt
1265 if (nstart_seq.eq.0) nstart_seq=nnt
1266 if(me.eq.king.or..not.out1file) &
1267 write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1268 ' nstart_seq=',nstart_seq !,"242343453254"
1270 !--- Zscore rms -------
1271 if (nz_start.eq.0) nz_start=nnt
1272 if (nz_end.eq.0 .and. nsup.gt.0) then
1274 else if (nz_end.eq.0) then
1277 if(me.eq.king.or..not.out1file)then
1278 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1279 write (iout,*) 'IZ_SC=',iz_sc
1281 !----------------------
1284 if (.not.pdbref) then
1285 call read_angles(inp,*38)
1287 38 write (iout,'(a)') 'Error reading reference structure.'
1289 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1290 stop 'Error reading reference structure'
1294 !zscore call geom_to_var(nvar,coord_exp_zs(1,1))
1301 cref(j,i,kkk)=c(j,i)
1304 call contact(.true.,ncont_ref,icont_ref,co)
1306 ! write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1308 if (constr_dist.gt.0) call read_dist_constr
1309 write (iout,*) "After read_dist_constr nhpb",nhpb
1310 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1312 if(me.eq.king.or..not.out1file) &
1313 write (iout,*) 'Contact order:',co
1315 if(me.eq.king.or..not.out1file) &
1316 write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1319 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1321 if(me.eq.king.or..not.out1file) &
1322 write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1323 icont_ref(1,i),' ',&
1324 restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1328 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1329 .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1330 modecalc.ne.10) then
1331 ! If input structure hasn't been supplied from the PDB file read or generate
1333 if (iranconf.eq.0 .and. .not. extconf) then
1334 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1335 write (iout,'(a)') 'Initial geometry will be read in.'
1337 read(inp,'(8f10.5)',end=36,err=36) &
1338 ((c(l,k),l=1,3),k=1,nres),&
1339 ((c(l,k+nres),l=1,3),k=nnt,nct)
1340 write (iout,*) "Exit READ_CART"
1341 write (iout,'(8f10.5)') &
1342 ((c(l,k),l=1,3),k=1,nres)
1343 write (iout,'(8f10.5)') &
1344 ((c(l,k+nres),l=1,3),k=nnt,nct)
1345 call int_from_cart1(.true.)
1346 write (iout,*) "Finish INT_TO_CART"
1349 dc(j,i)=c(j,i+1)-c(j,i)
1350 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1354 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1356 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1357 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1363 call read_angles(inp,*36)
1366 36 write (iout,'(a)') 'Error reading angle file.'
1368 call mpi_finalize( MPI_COMM_WORLD,IERR )
1370 stop 'Error reading angle file.'
1372 else if (extconf) then
1373 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1374 write (iout,'(a)') 'Extended chain initial geometry.'
1376 theta(i)=90d0*deg2rad
1379 phi(i)=180d0*deg2rad
1382 alph(i)=110d0*deg2rad
1385 omeg(i)=-120d0*deg2rad
1386 if (itype(i,1).le.0) omeg(i)=-omeg(i)
1389 if(me.eq.king.or..not.out1file) &
1390 write (iout,'(a)') 'Random-generated initial geometry.'
1394 if (me.eq.king .or. fg_rank.eq.0 .and. &
1395 ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1399 call gen_rand_conf(itmp,*30)
1401 30 write (iout,*) 'Failed to generate random conformation',&
1403 write (*,*) 'Processor:',me,&
1404 ' Failed to generate random conformation',&
1414 write (iout,'(a,i3,a)') 'Processor:',me,&
1415 ' error in generating random conformation.'
1416 write (*,'(a,i3,a)') 'Processor:',me,&
1417 ' error in generating random conformation.'
1420 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1426 call gen_rand_conf(itmp,*335)
1428 335 write (iout,*) 'Failed to generate random conformation',&
1430 write (*,*) 'Failed to generate random conformation',&
1433 write (iout,'(a,i3,a)') 'Processor:',me,&
1434 ' error in generating random conformation.'
1435 write (*,'(a,i3,a)') 'Processor:',me,&
1436 ' error in generating random conformation.'
1441 elseif (modecalc.eq.4) then
1442 read (inp,'(a)') intinname
1443 open (intin,file=intinname,status='old',err=333)
1444 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1445 write (iout,'(a)') 'intinname',intinname
1446 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1448 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1450 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1452 stop 'Error opening angle file.'
1456 ! Generate distance constraints, if the PDB structure is to be regularized.
1457 if (nthread.gt.0) then
1458 call read_threadbase
1461 if (me.eq.king .or. .not. out1file) &
1463 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1464 write (iout,'(/a,i3,a)') &
1465 'The chain contains',ns,' disulfide-bridging cysteines.'
1466 write (iout,'(20i4)') (iss(i),i=1,ns)
1468 write(iout,*)"Running with dynamic disulfide-bond formation"
1470 write (iout,'(/a/)') 'Pre-formed links are:'
1476 if (me.eq.king.or..not.out1file) &
1477 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1478 restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1484 if (ns.gt.0.and.dyn_ss) then
1488 forcon(i-nss)=forcon(i)
1495 dyn_ss_mask(iss(i))=.true.
1498 if (i2ndstr.gt.0) call secstrp2dihc
1499 ! call geom_to_var(nvar,x)
1500 ! call etotal(energia(0))
1501 ! call enerprint(energia(0))
1502 ! call briefout(0,etot)
1504 !d write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1505 !d write (iout,'(a)') 'Variable list:'
1506 !d write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1508 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1509 write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1510 'Processor',myrank,': end reading molecular data.'
1513 end subroutine molread
1514 !-----------------------------------------------------------------------------
1515 !-----------------------------------------------------------------------------