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=480) :: 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,'WBOND',wbond,1.0D0)
773 call reada(weightcard,'WTOR',wtor,1.0D0)
774 call reada(weightcard,'WTORD',wtor_d,1.0D0)
775 call reada(weightcard,'WSHIELD',wshield,0.05D0)
776 call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
777 call reada(weightcard,'WLT',wliptran,1.0D0)
778 call reada(weightcard,'WTUBE',wtube,1.0d0)
779 call reada(weightcard,'WANG',wang,1.0D0)
780 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
781 call reada(weightcard,'SCAL14',scal14,0.4D0)
782 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
783 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
784 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
785 call reada(weightcard,'TEMP0',temp0,300.0d0)
786 if (index(weightcard,'SOFT').gt.0) ipot=6
787 call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0)
788 ! 12/1/95 Added weight for the multi-body term WCORR
789 call reada(weightcard,'WCORRH',wcorr,1.0D0)
790 if (wcorr4.gt.0.0d0) wcorr=wcorr4
810 if(me.eq.king.or..not.out1file) &
811 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
812 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
814 10 format (/'Energy-term weights (unscaled):'// &
815 'WSCC= ',f10.6,' (SC-SC)'/ &
816 'WSCP= ',f10.6,' (SC-p)'/ &
817 'WELEC= ',f10.6,' (p-p electr)'/ &
818 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
819 'WBOND= ',f10.6,' (stretching)'/ &
820 'WANG= ',f10.6,' (bending)'/ &
821 'WSCLOC= ',f10.6,' (SC local)'/ &
822 'WTOR= ',f10.6,' (torsional)'/ &
823 'WTORD= ',f10.6,' (double torsional)'/ &
824 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
825 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
826 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
827 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
828 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
829 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
830 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
831 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
832 'WTURN6= ',f10.6,' (turns, 6th order)')
833 if(me.eq.king.or..not.out1file)then
834 if (wcorr4.gt.0.0d0) then
835 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
836 'between contact pairs of peptide groups'
837 write (iout,'(2(a,f5.3/))') &
838 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
839 'Range of quenching the correlation terms:',2*delt_corr
840 else if (wcorr.gt.0.0d0) then
841 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
842 'between contact pairs of peptide groups'
844 write (iout,'(a,f8.3)') &
845 'Scaling factor of 1,4 SC-p interactions:',scal14
846 write (iout,'(a,f8.3)') &
847 'General scaling factor of SC-p interactions:',scalscp
849 r0_corr=cutoff_corr-delt_corr
851 aad(i,1)=scalscp*aad(i,1)
852 aad(i,2)=scalscp*aad(i,2)
853 bad(i,1)=scalscp*bad(i,1)
854 bad(i,2)=scalscp*bad(i,2)
856 call rescale_weights(t_bath)
857 if(me.eq.king.or..not.out1file) &
858 write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
859 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
861 22 format (/'Energy-term weights (scaled):'// &
862 'WSCC= ',f10.6,' (SC-SC)'/ &
863 'WSCP= ',f10.6,' (SC-p)'/ &
864 'WELEC= ',f10.6,' (p-p electr)'/ &
865 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
866 'WBOND= ',f10.6,' (stretching)'/ &
867 'WANG= ',f10.6,' (bending)'/ &
868 'WSCLOC= ',f10.6,' (SC local)'/ &
869 'WTOR= ',f10.6,' (torsional)'/ &
870 'WTORD= ',f10.6,' (double torsional)'/ &
871 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
872 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
873 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
874 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
875 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
876 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
877 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
878 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
879 'WTURN6= ',f10.6,' (turns, 6th order)')
880 if(me.eq.king.or..not.out1file) &
881 write (iout,*) "Reference temperature for weights calculation:",&
883 call reada(weightcard,"D0CM",d0cm,3.78d0)
884 call reada(weightcard,"AKCM",akcm,15.1d0)
885 call reada(weightcard,"AKTH",akth,11.0d0)
886 call reada(weightcard,"AKCT",akct,12.0d0)
887 call reada(weightcard,"V1SS",v1ss,-1.08d0)
888 call reada(weightcard,"V2SS",v2ss,7.61d0)
889 call reada(weightcard,"V3SS",v3ss,13.7d0)
890 call reada(weightcard,"EBR",ebr,-5.50D0)
891 call reada(weightcard,"ATRISS",atriss,0.301D0)
892 call reada(weightcard,"BTRISS",btriss,0.021D0)
893 call reada(weightcard,"CTRISS",ctriss,1.001D0)
894 call reada(weightcard,"DTRISS",dtriss,1.001D0)
895 call reada(weightcard,"SSSCALE",ssscale,1.0D0)
896 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
898 call reada(weightcard,"HT",Ht,0.0D0)
900 ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
901 Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
902 akcm=akcm/wsc*ssscale
903 akth=akth/wsc*ssscale
904 akct=akct/wsc*ssscale
905 v1ss=v1ss/wsc*ssscale
906 v2ss=v2ss/wsc*ssscale
907 v3ss=v3ss/wsc*ssscale
909 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
912 if(me.eq.king.or..not.out1file) then
913 write (iout,*) "Parameters of the SS-bond potential:"
914 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
916 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
917 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
918 write (iout,*)" HT",Ht
919 print *,'indpdb=',indpdb,' pdbref=',pdbref
921 if (indpdb.gt.0 .or. pdbref) then
922 read(inp,'(a)') pdbfile
923 if(me.eq.king.or..not.out1file) &
924 write (iout,'(2a)') 'PDB data will be read from file ',&
925 pdbfile(:ilen(pdbfile))
926 open(ipdbin,file=pdbfile,status='old',err=33)
928 33 write (iout,'(a)') 'Error opening PDB file.'
931 ! print *,'Begin reading pdb data'
933 ! print *,'Finished reading pdb data'
934 if(me.eq.king.or..not.out1file) &
935 write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
936 ' nstart_sup=',nstart_sup !,"ergwergewrgae"
937 !el if(.not.allocated(itype_pdb))
938 allocate(itype_pdb(nres))
940 itype_pdb(i)=itype(i,1)
944 nct=nstart_sup+nsup-1
945 !el if(.not.allocated(icont_ref))
946 allocate(icont_ref(2,12*nres)) ! maxcont=12*maxres
947 call contact(.false.,ncont_ref,icont_ref,co)
950 if(me.eq.king.or..not.out1file) &
951 write(iout,*)'Adding sidechains'
955 if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
958 do while (fail.and.nsi.le.maxsi)
959 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
962 if(fail) write(iout,*)'Adding sidechain failed for res ',&
963 i,' after ',nsi,' trials'
969 if (indpdb.eq.0) then
971 allocate(sequence(maxres,5))
974 ! Read sequence if not taken from the pdb file.
976 read (inp,*) nres_molec(molec)
977 ! print *,'nres=',nres
978 if (iscode.gt.0) then
979 read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
981 read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
983 ! Convert sequence to numeric code
985 do i=1,nres_molec(molec)
986 itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
990 ! Read sequence if not taken from the pdb file.
992 read (inp,*) nres_molec(molec)
993 ! print *,'nres=',nres
994 ! allocate(sequence(maxres,5))
995 ! if (iscode.gt.0) then
996 read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
997 ! Convert sequence to numeric code
999 do i=1,nres_molec(molec)
1000 istype(i)=sugarcode(sequence(i,molec)(1:1),i)
1001 itype(i,1)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
1006 ! Read sequence if not taken from the pdb file.
1008 read (inp,*) nres_molec(molec)
1009 ! print *,'nres=',nres
1010 read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1011 ! Convert sequence to numeric code
1012 print *,nres_molec(molec)
1013 do i=1,nres_molec(molec)
1014 itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1019 nres=nres+nres_molec(i)
1020 print *,nres_molec(i)
1023 ! Assign initial virtual bond lengths
1024 !elwrite(iout,*) "test_alloc"
1025 if(.not.allocated(vbld)) allocate(vbld(2*nres))
1026 !elwrite(iout,*) "test_alloc"
1027 if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1028 !elwrite(iout,*) "test_alloc"
1034 vbld(i+nres)=dsc(iabs(itype(i,1)))
1035 vbld_inv(i+nres)=dsc_inv(iabs(itype(i,1)))
1036 ! write (iout,*) "i",i," itype",itype(i,1),
1037 ! & " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1041 ! print '(20i4)',(itype(i,1),i=1,nres)
1042 !----------------------------
1043 !el reallocate tables
1046 ! c_alloc(j,i)=c(j,i)
1047 ! dc_alloc(j,i)=dc(j,i)
1051 !elwrite(iout,*) "itype",i,itype(i,1)
1052 ! itype_alloc(i)=itype(i,1)
1058 ! allocate(c(3,2*nres+4))
1059 ! allocate(dc(3,0:2*nres+2))
1060 ! allocate(itype(nres+2))
1061 allocate(itel(nres+2))
1066 ! c(j,i)=c_alloc(j,i)
1067 ! dc(j,i)=dc_alloc(j,i)
1071 ! itype(i,1)=itype_alloc(i)
1074 !--------------------------
1077 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1079 if (itype(i,1).eq.ntyp1) then
1083 else if (iabs(itype(i+1,1)).ne.20) then
1085 else if (iabs(itype(i,1)).ne.20) then
1092 if(me.eq.king.or..not.out1file)then
1093 write (iout,*) "ITEL"
1095 write (iout,*) i,itype(i,1),itel(i)
1097 print *,'Call Read_Bridge.'
1100 !--------------------------------
1101 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1102 call alloc_geo_arrays
1103 call alloc_ener_arrays
1104 !--------------------------------
1105 ! 8/13/98 Set limits to generating the dihedral angles
1110 read (inp,*) ndih_constr
1111 if (ndih_constr.gt.0) then
1112 allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1113 allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1116 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1117 if(me.eq.king.or..not.out1file)then
1119 'There are',ndih_constr,' constraints on phi angles.'
1121 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1125 phi0(i)=deg2rad*phi0(i)
1126 drange(i)=deg2rad*drange(i)
1128 if(me.eq.king.or..not.out1file) &
1129 write (iout,*) 'FTORS',ftors
1132 phibound(1,ii) = phi0(i)-drange(i)
1133 phibound(2,ii) = phi0(i)+drange(i)
1136 if (with_theta_constr) then
1137 !C with_theta_constr is keyword allowing for occurance of theta constrains
1138 read (inp,*) ntheta_constr
1139 !C ntheta_constr is the number of theta constrains
1140 if (ntheta_constr.gt.0) then
1141 !C read (inp,*) ftors
1142 allocate(itheta_constr(ntheta_constr))
1143 allocate(theta_constr0(ntheta_constr))
1144 allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1145 read (inp,*) (itheta_constr(i),theta_constr0(i), &
1146 theta_drange(i),for_thet_constr(i), &
1148 !C the above code reads from 1 to ntheta_constr
1149 !C itheta_constr(i) residue i for which is theta_constr
1150 !C theta_constr0 the global minimum value
1151 !C theta_drange is range for which there is no energy penalty
1152 !C for_thet_constr is the force constant for quartic energy penalty
1154 if(me.eq.king.or..not.out1file)then
1156 'There are',ntheta_constr,' constraints on phi angles.'
1157 do i=1,ntheta_constr
1158 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1163 do i=1,ntheta_constr
1164 theta_constr0(i)=deg2rad*theta_constr0(i)
1165 theta_drange(i)=deg2rad*theta_drange(i)
1167 !C if(me.eq.king.or..not.out1file)
1168 !C & write (iout,*) 'FTORS',ftors
1169 !C do i=1,ntheta_constr
1170 !C ii = itheta_constr(i)
1171 !C thetabound(1,ii) = phi0(i)-drange(i)
1172 !C thetabound(2,ii) = phi0(i)+drange(i)
1174 endif ! ntheta_constr.gt.0
1175 endif! with_theta_constr
1179 if (me.eq.king) then
1181 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1183 write (iout,'(a3,i5,2f10.1)') &
1184 restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1190 !d print *,'NNT=',NNT,' NCT=',NCT
1191 if (itype(1,1).eq.ntyp1) nnt=2
1192 if (itype(nres,1).eq.ntyp1) nct=nct-1
1194 if(me.eq.king.or..not.out1file) &
1195 write (iout,'(a,i3)') 'nsup=',nsup
1197 if (nsup.le.(nct-nnt+1)) then
1198 do i=0,nct-nnt+1-nsup
1199 if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1204 write (iout,'(a)') &
1205 'Error - sequences to be superposed do not match.'
1208 do i=0,nsup-(nct-nnt+1)
1209 if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1211 nstart_sup=nstart_sup+i
1216 write (iout,'(a)') &
1217 'Error - sequences to be superposed do not match.'
1220 if (nsup.eq.0) nsup=nct-nnt
1221 if (nstart_sup.eq.0) nstart_sup=nnt
1222 if (nstart_seq.eq.0) nstart_seq=nnt
1223 if(me.eq.king.or..not.out1file) &
1224 write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1225 ' nstart_seq=',nstart_seq !,"242343453254"
1227 !--- Zscore rms -------
1228 if (nz_start.eq.0) nz_start=nnt
1229 if (nz_end.eq.0 .and. nsup.gt.0) then
1231 else if (nz_end.eq.0) then
1234 if(me.eq.king.or..not.out1file)then
1235 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1236 write (iout,*) 'IZ_SC=',iz_sc
1238 !----------------------
1241 if (.not.pdbref) then
1242 call read_angles(inp,*38)
1244 38 write (iout,'(a)') 'Error reading reference structure.'
1246 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1247 stop 'Error reading reference structure'
1251 !zscore call geom_to_var(nvar,coord_exp_zs(1,1))
1258 cref(j,i,kkk)=c(j,i)
1261 call contact(.true.,ncont_ref,icont_ref,co)
1263 ! write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1265 if (constr_dist.gt.0) call read_dist_constr
1266 write (iout,*) "After read_dist_constr nhpb",nhpb
1267 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1269 if(me.eq.king.or..not.out1file) &
1270 write (iout,*) 'Contact order:',co
1272 if(me.eq.king.or..not.out1file) &
1273 write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1276 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1278 if(me.eq.king.or..not.out1file) &
1279 write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1280 icont_ref(1,i),' ',&
1281 restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1285 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1286 .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1287 modecalc.ne.10) then
1288 ! If input structure hasn't been supplied from the PDB file read or generate
1290 if (iranconf.eq.0 .and. .not. extconf) then
1291 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1292 write (iout,'(a)') 'Initial geometry will be read in.'
1294 read(inp,'(8f10.5)',end=36,err=36) &
1295 ((c(l,k),l=1,3),k=1,nres),&
1296 ((c(l,k+nres),l=1,3),k=nnt,nct)
1297 write (iout,*) "Exit READ_CART"
1298 write (iout,'(8f10.5)') &
1299 ((c(l,k),l=1,3),k=1,nres),&
1300 ((c(l,k+nres),l=1,3),k=nnt,nct)
1301 call int_from_cart1(.true.)
1302 write (iout,*) "Finish INT_TO_CART"
1305 dc(j,i)=c(j,i+1)-c(j,i)
1306 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1310 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1312 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1313 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1319 call read_angles(inp,*36)
1322 36 write (iout,'(a)') 'Error reading angle file.'
1324 call mpi_finalize( MPI_COMM_WORLD,IERR )
1326 stop 'Error reading angle file.'
1328 else if (extconf) then
1329 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1330 write (iout,'(a)') 'Extended chain initial geometry.'
1332 theta(i)=90d0*deg2rad
1335 phi(i)=180d0*deg2rad
1338 alph(i)=110d0*deg2rad
1340 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1342 omeg(i)=-120d0*deg2rad
1343 if (itype(i,1).le.0) omeg(i)=-omeg(i)
1346 if(me.eq.king.or..not.out1file) &
1347 write (iout,'(a)') 'Random-generated initial geometry.'
1351 if (me.eq.king .or. fg_rank.eq.0 .and. &
1352 ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1356 call gen_rand_conf(itmp,*30)
1358 30 write (iout,*) 'Failed to generate random conformation',&
1360 write (*,*) 'Processor:',me,&
1361 ' Failed to generate random conformation',&
1371 write (iout,'(a,i3,a)') 'Processor:',me,&
1372 ' error in generating random conformation.'
1373 write (*,'(a,i3,a)') 'Processor:',me,&
1374 ' error in generating random conformation.'
1377 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1383 call gen_rand_conf(itmp,*335)
1385 335 write (iout,*) 'Failed to generate random conformation',&
1387 write (*,*) 'Failed to generate random conformation',&
1390 write (iout,'(a,i3,a)') 'Processor:',me,&
1391 ' error in generating random conformation.'
1392 write (*,'(a,i3,a)') 'Processor:',me,&
1393 ' error in generating random conformation.'
1398 elseif (modecalc.eq.4) then
1399 read (inp,'(a)') intinname
1400 open (intin,file=intinname,status='old',err=333)
1401 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1402 write (iout,'(a)') 'intinname',intinname
1403 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1405 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1407 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1409 stop 'Error opening angle file.'
1413 ! Generate distance constraints, if the PDB structure is to be regularized.
1414 if (nthread.gt.0) then
1415 call read_threadbase
1418 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1419 if (me.eq.king .or. .not. out1file) &
1421 !elwrite (iout,*)"alph(i)*rad2deg",(alph(i)*rad2deg, i=1,nres)
1422 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1423 write (iout,'(/a,i3,a)') &
1424 'The chain contains',ns,' disulfide-bridging cysteines.'
1425 write (iout,'(20i4)') (iss(i),i=1,ns)
1427 write(iout,*)"Running with dynamic disulfide-bond formation"
1429 write (iout,'(/a/)') 'Pre-formed links are:'
1435 if (me.eq.king.or..not.out1file) &
1436 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1437 restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1443 if (ns.gt.0.and.dyn_ss) then
1447 forcon(i-nss)=forcon(i)
1454 dyn_ss_mask(iss(i))=.true.
1457 if (i2ndstr.gt.0) call secstrp2dihc
1458 ! call geom_to_var(nvar,x)
1459 ! call etotal(energia(0))
1460 ! call enerprint(energia(0))
1461 ! call briefout(0,etot)
1463 !d write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1464 !d write (iout,'(a)') 'Variable list:'
1465 !d write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1467 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1468 write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1469 'Processor',myrank,': end reading molecular data.'
1472 end subroutine molread
1473 !-----------------------------------------------------------------------------
1474 !-----------------------------------------------------------------------------