12 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 !-----------------------------------------------------------------------------
19 !-----------------------------------------------------------------------------
25 ! implicit real*8 (a-h,o-z)
26 ! include 'DIMENSIONS'
27 ! include 'DIMENSIONS.ZSCOPT'
31 ! include 'COMMON.MPI'
33 character(len=3) :: liczba
35 ! include 'COMMON.IOUNITS'
36 integer :: lenpre,lenpot !,ilen
42 call mygetenv('PREFIX',prefix)
43 call mygetenv('SCRATCHDIR',scratchdir)
44 call mygetenv('POT',pot)
47 call mygetenv('POT',pot)
48 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
49 ! Get the names and open the input files
50 open (1,file=prefix(:ilen(prefix))//'.inp',status='old')
51 ! Get parameter filenames and open the parameter files.
52 call mygetenv('BONDPAR',bondname)
53 open (ibond,file=bondname,status='old')
54 call mygetenv('THETPAR',thetname)
55 open (ithep,file=thetname,status='old')
56 call mygetenv('ROTPAR',rotname)
57 open (irotam,file=rotname,status='old')
58 call mygetenv('TORPAR',torname)
59 open (itorp,file=torname,status='old')
60 call mygetenv('TORDPAR',tordname)
61 open (itordp,file=tordname,status='old')
62 call mygetenv('FOURIER',fouriername)
63 open (ifourier,file=fouriername,status='old')
64 call mygetenv('SCCORPAR',sccorname)
65 open (isccor,file=sccorname,status='old')
66 call mygetenv('ELEPAR',elename)
67 open (ielep,file=elename,status='old')
68 call mygetenv('SIDEPAR',sidename)
69 open (isidep,file=sidename,status='old')
70 call mygetenv('SIDEP',sidepname)
71 open (isidep1,file=sidepname,status="old")
74 ! 8/9/01 In the newest version SCp interaction constants are read from a file
75 ! Use -DOLDSCP to use hard-coded constants instead.
77 call mygetenv('SCPPAR',scpname)
78 open (iscpp,file=scpname,status='old')
81 if (MyID.eq.BossID) then
85 print *,'OpenUnits: processor',MyRank
86 call numstr(MyRank,liczba)
87 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//liczba
89 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
91 open(iout,file=outname,status='unknown')
92 write (iout,'(80(1h-))')
93 write (iout,'(30x,a)') "FILE ASSIGNMENT"
94 write (iout,'(80(1h-))')
95 write (iout,*) "Input file : ",&
96 prefix(:ilen(prefix))//'.inp'
97 write (iout,*) "Output file : ",&
98 outname(:ilen(outname))
100 write (iout,*) "Sidechain potential file : ",&
101 sidename(:ilen(sidename))
103 write (iout,*) "SCp potential file : ",&
104 scpname(:ilen(scpname))
106 write (iout,*) "Electrostatic potential file : ",&
107 elename(:ilen(elename))
108 write (iout,*) "Cumulant coefficient file : ",&
109 fouriername(:ilen(fouriername))
110 write (iout,*) "Torsional parameter file : ",&
111 torname(:ilen(torname))
112 write (iout,*) "Double torsional parameter file : ",&
113 tordname(:ilen(tordname))
114 write (iout,*) "Backbone-rotamer parameter file : ",&
115 sccorname(:ilen(sccorname))
116 write (iout,*) "Bond & inertia constant file : ",&
117 bondname(:ilen(bondname))
118 write (iout,*) "Bending parameter file : ",&
119 thetname(:ilen(thetname))
120 write (iout,*) "Rotamer parameter file : ",&
121 rotname(:ilen(rotname))
122 write (iout,'(80(1h-))')
125 end subroutine openunits
126 !-----------------------------------------------------------------------------
128 !-----------------------------------------------------------------------------
129 subroutine molread(*)
131 ! Read molecular data.
134 use geometry_data, only:nres,deg2rad,c,dc
135 use control_data, only:iscode
136 use control, only:rescode,setup_var,init_int_table
137 use geometry, only:alloc_geo_arrays
138 use energy, only:alloc_ener_arrays
139 ! implicit real*8 (a-h,o-z)
140 ! include 'DIMENSIONS'
141 ! include 'DIMENSIONS.ZSCOPT'
142 ! include 'COMMON.IOUNITS'
143 ! include 'COMMON.GEO'
144 ! include 'COMMON.VAR'
145 ! include 'COMMON.INTERACT'
146 ! include 'COMMON.LOCAL'
147 ! include 'COMMON.NAMES'
148 ! include 'COMMON.CHAIN'
149 ! include 'COMMON.FFIELD'
150 ! include 'COMMON.SBRIDGE'
151 ! include 'COMMON.TORCNSTR'
152 ! include 'COMMON.CONTROL'
153 character(len=4),dimension(:),allocatable :: sequence !(nres)
154 !el integer :: rescode
155 !el real(kind=8) :: x(maxvar)
156 character(len=320) :: controlcard !,ucase
157 integer,dimension(nres) :: itype_pdb !(maxres)
158 integer :: i,j,i1,i2,it1,it2
159 real(kind=8) :: scalscp
160 !el logical :: seq_comp
161 call card_concat(controlcard,.true.)
162 call reada(controlcard,'SCAL14',scal14,0.4d0)
163 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
164 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
165 call reada(controlcard,'TEMP0',temp0,300.0d0) !el
166 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
167 r0_corr=cutoff_corr-delt_corr
168 call readi(controlcard,"NRES",nres,0)
169 allocate(sequence(nres+1))
170 !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
171 call alloc_geo_arrays
172 call alloc_ener_arrays
173 ! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
174 !----------------------------
175 allocate(c(3,2*nres+2))
176 allocate(dc(3,0:2*nres+2))
177 allocate(itype(nres+2))
178 allocate(itel(nres+2))
191 !--------------------------
193 iscode=index(controlcard,"ONE_LETTER")
195 write (iout,*) "Error: no residues in molecule"
198 if (nres.gt.maxres) then
199 write (iout,*) "Error: too many residues",nres,maxres
201 write(iout,*) 'nres=',nres
202 ! Read sequence of the protein
203 if (iscode.gt.0) then
204 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
206 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
208 ! Convert sequence to numeric code
210 itype(i)=rescode(i,sequence(i),iscode)
212 write (iout,*) "Numeric code:"
213 write (iout,'(20i4)') (itype(i),i=1,nres)
216 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
218 if (itype(i).eq.ntyp1) then
222 else if (iabs(itype(i+1)).ne.20) then
224 else if (iabs(itype(i)).ne.20) then
231 write (iout,*) "ITEL"
233 write (iout,*) i,itype(i),itel(i)
237 if (with_dihed_constr) then
239 read (inp,*) ndih_constr
240 if (ndih_constr.gt.0) then
242 write (iout,*) 'FTORS',ftors
243 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
245 'There are',ndih_constr,' constraints on phi angles.'
247 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
250 phi0(i)=deg2rad*phi0(i)
251 drange(i)=deg2rad*drange(i)
259 if (itype(1).eq.ntyp1) nnt=2
260 if (itype(nres).eq.ntyp1) nct=nct-1
261 write(iout,*) 'NNT=',NNT,' NCT=',NCT
265 write (iout,'(/a,i3,a)') 'The chain contains',ns,&
266 ' disulfide-bridging cysteines.'
267 write (iout,'(20i4)') (iss(i),i=1,ns)
268 write (iout,'(/a/)') 'Pre-formed links are:'
274 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
275 restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',&
276 dhpb(i),ebr,forcon(i)
281 end subroutine molread
282 !-----------------------------------------------------------------------------
284 !-----------------------------------------------------------------------------
285 subroutine parmread(iparm,*)
290 ! Read the parameters of the probability distributions of the virtual-bond
291 ! valence angles and the side chains and energy parameters.
297 use control_data, only: maxtor,maxterm,maxlor,maxterm_sccor,&
298 maxtermd_1,maxtermd_2,maxthetyp,maxthetyp1
302 use io_config, only: printmat
303 use control, only: getenv_loc
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'DIMENSIONS.ZSCOPT'
313 ! include 'DIMENSIONS.FREE'
314 ! include 'COMMON.IOUNITS'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.INTERACT'
317 ! include 'COMMON.GEO'
318 ! include 'COMMON.LOCAL'
319 ! include 'COMMON.TORSION'
320 ! include 'COMMON.FFIELD'
321 ! include 'COMMON.NAMES'
322 ! include 'COMMON.SBRIDGE'
323 ! include 'COMMON.WEIGHTS'
324 ! include 'COMMON.ENEPS'
325 ! include 'COMMON.SCCOR'
326 ! include 'COMMON.SCROT'
327 ! include 'COMMON.FREE'
328 character(len=1) :: t1,t2,t3
329 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
330 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
332 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
333 character(len=800) :: controlcard
334 character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
335 tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
339 character(len=16) :: key
341 !el real(kind=8) :: ip,mp
342 real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,&
343 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
344 real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
346 integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n
347 integer :: nlobi,iblock,maxinter,iscprol
351 ! Set LPRINT=.TRUE. for debugging
352 dwa16=2.0d0**(1.0d0/6.0d0)
355 ! Assign virtual-bond length
358 vblinv2=vblinv*vblinv
360 call card_concat(controlcard,.true.)
363 allocate(ww(max_eneW))
365 key = wname(i)(:ilen(wname(i)))
366 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
369 write (iout,*) "iparm",iparm," myparm",myparm
370 ! If reading not own parameters, skip assignment
372 if (iparm.eq.myparm .or. .not.separate_parset) then
375 ! Setup weights for UNRES
398 allocate(weights(n_ene))
413 weights(15)=0 !wstrain !
414 weights(16)=0 !wvdwpp !
416 weights(18)=0 !scal14 !
419 call card_concat(controlcard,.false.)
421 ! Return if not own parameters
423 if (iparm.ne.myparm .and. separate_parset) return
425 call reads(controlcard,"BONDPAR",bondname_t,bondname)
426 open (ibond,file=bondname_t,status='old')
428 call reads(controlcard,"THETPAR",thetname_t,thetname)
429 open (ithep,file=thetname_t,status='old')
431 call reads(controlcard,"ROTPAR",rotname_t,rotname)
432 open (irotam,file=rotname_t,status='old')
434 call reads(controlcard,"TORPAR",torname_t,torname)
435 open (itorp,file=torname_t,status='old')
437 call reads(controlcard,"TORDPAR",tordname_t,tordname)
438 open (itordp,file=tordname_t,status='old')
440 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
441 open (isccor,file=sccorname_t,status='old')
443 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
444 open (ifourier,file=fouriername_t,status='old')
446 call reads(controlcard,"ELEPAR",elename_t,elename)
447 open (ielep,file=elename_t,status='old')
449 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
450 open (isidep,file=sidename_t,status='old')
452 call reads(controlcard,"SCPPAR",scpname_t,scpname)
453 open (iscpp,file=scpname_t,status='old')
455 write (iout,*) "Parameter set:",iparm
456 write (iout,*) "Energy-term weights:"
458 write (iout,'(a16,f10.5)') wname(i),ww(i)
460 write (iout,*) "Sidechain potential file : ",&
461 sidename_t(:ilen(sidename_t))
463 write (iout,*) "SCp potential file : ",&
464 scpname_t(:ilen(scpname_t))
466 write (iout,*) "Electrostatic potential file : ",&
467 elename_t(:ilen(elename_t))
468 write (iout,*) "Cumulant coefficient file : ",&
469 fouriername_t(:ilen(fouriername_t))
470 write (iout,*) "Torsional parameter file : ",&
471 torname_t(:ilen(torname_t))
472 write (iout,*) "Double torsional parameter file : ",&
473 tordname_t(:ilen(tordname_t))
474 write (iout,*) "Backbone-rotamer parameter file : ",&
475 sccorname(:ilen(sccorname))
476 write (iout,*) "Bond & inertia constant file : ",&
477 bondname_t(:ilen(bondname_t))
478 write (iout,*) "Bending parameter file : ",&
479 thetname_t(:ilen(thetname_t))
480 write (iout,*) "Rotamer parameter file : ",&
481 rotname_t(:ilen(rotname_t))
484 ! Read the virtual-bond parameters, masses, and moments of inertia
485 ! and Stokes' radii of the peptide group and side chains
487 allocate(dsc(ntyp1)) !(ntyp1)
488 allocate(dsc_inv(ntyp1)) !(ntyp1)
489 allocate(nbondterm(ntyp)) !(ntyp)
490 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
491 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
492 !el allocate(msc(ntyp+1)) !(ntyp+1)
493 !el allocate(isc(ntyp+1)) !(ntyp+1)
494 !el allocate(restok(ntyp+1)) !(ntyp+1)
495 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
498 read (ibond,*) vbldp0,akp
501 read (ibond,*) vbldsc0(1,i),aksc(1,i)
502 dsc(i) = vbldsc0(1,i)
506 dsc_inv(i)=1.0D0/dsc(i)
510 read (ibond,*) ijunk,vbldp0,akp,rjunk
512 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
514 dsc(i) = vbldsc0(1,i)
518 dsc_inv(i)=1.0D0/dsc(i)
523 write(iout,'(/a/)')"Force constants virtual bonds:"
524 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
526 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
528 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
529 vbldsc0(1,i),aksc(1,i),abond0(1,i)
531 write (iout,'(13x,3f10.5)') &
532 vbldsc0(j,i),aksc(j,i),abond0(j,i)
536 !----------------------------------------------------
537 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
538 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
539 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
540 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
541 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
542 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
548 athet(j,i,ichir1,ichir2)=0.0D0
549 bthet(j,i,ichir1,ichir2)=0.0D0
563 !elwrite(iout,*) "parmread kontrol"
567 ! Read the parameters of the probability distribution/energy expression
568 ! of the virtual-bond valence angles theta
571 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
572 (bthet(j,i,1,1),j=1,2)
573 read (ithep,*) (polthet(j,i),j=0,3)
574 !elwrite(iout,*) "parmread kontrol in cryst_theta"
575 read (ithep,*) (gthet(j,i),j=1,3)
576 !elwrite(iout,*) "parmread kontrol in cryst_theta"
577 read (ithep,*) theta0(i),sig0(i),sigc0(i)
579 !elwrite(iout,*) "parmread kontrol in cryst_theta"
581 !elwrite(iout,*) "parmread kontrol in cryst_theta"
583 athet(1,i,1,-1)=athet(1,i,1,1)
584 athet(2,i,1,-1)=athet(2,i,1,1)
585 bthet(1,i,1,-1)=-bthet(1,i,1,1)
586 bthet(2,i,1,-1)=-bthet(2,i,1,1)
587 athet(1,i,-1,1)=-athet(1,i,1,1)
588 athet(2,i,-1,1)=-athet(2,i,1,1)
589 bthet(1,i,-1,1)=bthet(1,i,1,1)
590 bthet(2,i,-1,1)=bthet(2,i,1,1)
592 !elwrite(iout,*) "parmread kontrol in cryst_theta"
595 athet(1,i,-1,-1)=athet(1,-i,1,1)
596 athet(2,i,-1,-1)=-athet(2,-i,1,1)
597 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
598 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
599 athet(1,i,-1,1)=athet(1,-i,1,1)
600 athet(2,i,-1,1)=-athet(2,-i,1,1)
601 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
602 bthet(2,i,-1,1)=bthet(2,-i,1,1)
603 athet(1,i,1,-1)=-athet(1,-i,1,1)
604 athet(2,i,1,-1)=athet(2,-i,1,1)
605 bthet(1,i,1,-1)=bthet(1,-i,1,1)
606 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
611 polthet(j,i)=polthet(j,-i)
614 gthet(j,i)=gthet(j,-i)
617 !elwrite(iout,*) "parmread kontrol in cryst_theta"
619 !elwrite(iout,*) "parmread kontrol in cryst_theta"
622 ! & 'Parameters of the virtual-bond valence angles:'
623 ! write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
624 ! & ' ATHETA0 ',' A1 ',' A2 ',
627 ! write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
628 ! & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
630 ! write (iout,'(/a/9x,5a/79(1h-))')
631 ! & 'Parameters of the expression for sigma(theta_c):',
632 ! & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
633 ! & ' ALPH3 ',' SIGMA0C '
635 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
636 ! & (polthet(j,i),j=0,3),sigc0(i)
638 ! write (iout,'(/a/9x,5a/79(1h-))')
639 ! & 'Parameters of the second gaussian:',
640 ! & ' THETA0 ',' SIGMA0 ',' G1 ',
643 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
644 ! & sig0(i),(gthet(j,i),j=1,3)
647 'Parameters of the virtual-bond valence angles:'
648 write (iout,'(/a/9x,5a/79(1h-))') &
649 'Coefficients of expansion',&
650 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
651 ' b1*10^1 ',' b2*10^1 '
653 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
654 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
655 (10*bthet(j,i,1,1),j=1,2)
657 write (iout,'(/a/9x,5a/79(1h-))') &
658 'Parameters of the expression for sigma(theta_c):',&
659 ' alpha0 ',' alph1 ',' alph2 ',&
660 ' alhp3 ',' sigma0c '
662 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
663 (polthet(j,i),j=0,3),sigc0(i)
665 write (iout,'(/a/9x,5a/79(1h-))') &
666 'Parameters of the second gaussian:',&
667 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
670 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
671 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
676 ! Read the parameters of Utheta determined from ab initio surfaces
677 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
679 ! write (iout,*) "tu dochodze"
680 read (ithep,*) nthetyp,ntheterm,ntheterm2,&
681 ntheterm3,nsingle,ndouble
682 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
684 !----------------------------------------------------
685 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
686 allocate(aa0thet(-maxthetyp1:maxthetyp1,&
687 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
688 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
689 allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,&
690 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
691 !(maxtheterm,-maxthetyp1:maxthetyp1,&
692 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
693 allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
694 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
695 allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
696 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
697 allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
698 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
699 allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
700 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
701 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
702 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
703 allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
704 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
705 allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
706 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
707 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
708 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
711 read (ithep,*) (ithetyp(i),i=1,ntyp1)
713 ithetyp(i)=-ithetyp(-i)
715 ! write (iout,*) "tu dochodze"
717 do i=-maxthetyp,maxthetyp
718 do j=-maxthetyp,maxthetyp
719 do k=-maxthetyp,maxthetyp
720 aa0thet(i,j,k,iblock)=0.0d0
722 aathet(l,i,j,k,iblock)=0.0d0
726 bbthet(m,l,i,j,k,iblock)=0.0d0
727 ccthet(m,l,i,j,k,iblock)=0.0d0
728 ddthet(m,l,i,j,k,iblock)=0.0d0
729 eethet(m,l,i,j,k,iblock)=0.0d0
735 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
736 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
746 do j=-nthetyp,nthetyp
747 do k=-nthetyp,nthetyp
748 read (ithep,'(6a)') res1
749 read (ithep,*) aa0thet(i,j,k,iblock)
750 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
752 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
753 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
754 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
755 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
758 (((ffthet(llll,lll,ll,i,j,k,iblock),&
759 ffthet(lll,llll,ll,i,j,k,iblock),&
760 ggthet(llll,lll,ll,i,j,k,iblock),&
761 ggthet(lll,llll,ll,i,j,k,iblock),&
762 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
767 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
768 ! coefficients of theta-and-gamma-dependent terms are zero.
773 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
774 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
776 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
777 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
780 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
782 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
785 ! Substitution for D aminoacids from symmetry.
788 do j=-nthetyp,nthetyp
789 do k=-nthetyp,nthetyp
790 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
792 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
796 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
797 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
798 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
799 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
805 ffthet(llll,lll,ll,i,j,k,iblock)= &
806 ffthet(llll,lll,ll,-i,-j,-k,iblock)
807 ffthet(lll,llll,ll,i,j,k,iblock)= &
808 ffthet(lll,llll,ll,-i,-j,-k,iblock)
809 ggthet(llll,lll,ll,i,j,k,iblock)= &
810 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
811 ggthet(lll,llll,ll,i,j,k,iblock)= &
812 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
822 ! Control printout of the coefficients of virtual-bond-angle potentials
826 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
830 write (iout,'(//4a)') &
831 'Type ',onelett(i),onelett(j),onelett(k)
832 write (iout,'(//a,10x,a)') " l","a[l]"
833 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
834 write (iout,'(i2,1pe15.5)') &
835 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
837 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
838 "b",l,"c",l,"d",l,"e",l
840 write (iout,'(i2,4(1pe15.5))') m,&
841 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
842 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
846 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
847 "f+",l,"f-",l,"g+",l,"g-",l
850 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
851 ffthet(n,m,l,i,j,k,iblock),&
852 ffthet(m,n,l,i,j,k,iblock),&
853 ggthet(n,m,l,i,j,k,iblock),&
854 ggthet(m,n,l,i,j,k,iblock)
865 !-------------------------------------------
866 allocate(nlob(ntyp1)) !(ntyp1)
867 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
868 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
869 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
887 gaussc(l,k,j,i)=0.0D0
895 ! Read the parameters of the probability distribution/energy expression
896 ! of the side chains.
899 !c write (iout,*) "tu dochodze",i
900 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
904 dsc_inv(i)=1.0D0/dsc(i)
915 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
916 censc(1,1,-i)=censc(1,1,i)
917 censc(2,1,-i)=censc(2,1,i)
918 censc(3,1,-i)=-censc(3,1,i)
920 read (irotam,*) bsc(j,i)
921 read (irotam,*) (censc(k,j,i),k=1,3),&
922 ((blower(k,l,j),l=1,k),k=1,3)
923 censc(1,j,-i)=censc(1,j,i)
924 censc(2,j,-i)=censc(2,j,i)
925 censc(3,j,-i)=-censc(3,j,i)
926 ! BSC is amplitude of Gaussian
933 akl=akl+blower(k,m,j)*blower(l,m,j)
937 if (((k.eq.3).and.(l.ne.3)) &
938 .or.((l.eq.3).and.(k.ne.3))) then
939 gaussc(k,l,j,-i)=-akl
940 gaussc(l,k,j,-i)=-akl
952 write (iout,'(/a)') 'Parameters of side-chain local geometry'
956 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
957 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
958 ! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
959 ! write (iout,'(a,f10.4,4(16x,f10.4))')
960 ! & 'Center ',(bsc(j,i),j=1,nlobi)
961 ! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
962 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
963 'log h',(bsc(j,i),j=1,nlobi)
964 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
965 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
972 ! blower(k,l,j)=gaussc(ind,j,i)
977 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
978 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
985 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
986 ! added by Urszula Kozlowska 07/11/2007
988 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
996 read(irotam,*) sc_parmin(j,i)
1004 ! Read torsional parameters in old format
1006 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1008 read (itorp,*) ntortyp,nterm_old
1009 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1010 read (itorp,*) (itortyp(i),i=1,ntyp)
1012 !el from energy module--------
1013 allocate(v1(nterm_old,ntortyp,ntortyp))
1014 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1015 !el---------------------------
1021 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1027 write (iout,'(/a/)') 'Torsional constants:'
1030 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1031 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1039 ! Read torsional parameters
1041 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1043 read (itorp,*) ntortyp
1044 read (itorp,*) (itortyp(i),i=1,ntyp)
1045 write (iout,*) 'ntortyp',ntortyp
1047 !el from energy module---------
1048 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1049 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1051 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1052 allocate(vlor2(maxlor,ntortyp,ntortyp))
1053 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1054 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1056 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1057 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1058 !el---------------------------
1060 do i=-ntortyp,ntortyp
1061 do j=-ntortyp,ntortyp
1067 !el---------------------------
1071 itortyp(i)=-itortyp(-i)
1073 ! write (iout,*) 'ntortyp',ntortyp
1075 do j=-ntortyp+1,ntortyp-1
1076 read (itorp,*) nterm(i,j,iblock),&
1078 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1079 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1082 do k=1,nterm(i,j,iblock)
1083 read (itorp,*) kk,v1(k,i,j,iblock),&
1085 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1086 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1087 v0ij=v0ij+si*v1(k,i,j,iblock)
1090 do k=1,nlor(i,j,iblock)
1091 read (itorp,*) kk,vlor1(k,i,j),&
1092 vlor2(k,i,j),vlor3(k,i,j)
1093 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1096 v0(-i,-j,iblock)=v0ij
1103 write (iout,'(/a/)') 'Torsional constants:'
1106 write (iout,*) 'ityp',i,' jtyp',j
1107 write (iout,*) 'Fourier constants'
1108 do k=1,nterm(i,j,iblock)
1109 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1112 write (iout,*) 'Lorenz constants'
1113 do k=1,nlor(i,j,iblock)
1114 write (iout,'(3(1pe15.5))') &
1115 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1122 ! 6/23/01 Read parameters for double torsionals
1124 !el from energy module------------
1125 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1126 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1127 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1128 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1129 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1130 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1131 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1132 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1133 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1134 !---------------------------------
1138 do j=-ntortyp+1,ntortyp-1
1139 do k=-ntortyp+1,ntortyp-1
1140 read (itordp,'(3a1)') t1,t2,t3
1141 ! write (iout,*) "OK onelett",
1144 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1145 .or. t3.ne.toronelet(k)) then
1146 write (iout,*) "Error in double torsional parameter file",&
1149 call MPI_Finalize(Ierror)
1151 stop "Error in double torsional parameter file"
1153 read (itordp,*) ntermd_1(i,j,k,iblock),&
1154 ntermd_2(i,j,k,iblock)
1155 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1156 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1157 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1158 ntermd_1(i,j,k,iblock))
1159 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1160 ntermd_1(i,j,k,iblock))
1161 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1162 ntermd_1(i,j,k,iblock))
1163 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1164 ntermd_1(i,j,k,iblock))
1165 ! Martix of D parameters for one dimesional foureir series
1166 do l=1,ntermd_1(i,j,k,iblock)
1167 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1168 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1169 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1170 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1171 ! write(iout,*) "whcodze" ,
1172 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1174 read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1175 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1176 v2s(m,l,i,j,k,iblock),&
1177 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1178 ! Martix of D parameters for two dimesional fourier series
1179 do l=1,ntermd_2(i,j,k,iblock)
1181 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1182 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1183 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1184 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1193 write (iout,*) 'Constants for double torsionals'
1196 do j=-ntortyp+1,ntortyp-1
1197 do k=-ntortyp+1,ntortyp-1
1198 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1199 ' nsingle',ntermd_1(i,j,k,iblock),&
1200 ' ndouble',ntermd_2(i,j,k,iblock)
1202 write (iout,*) 'Single angles:'
1203 do l=1,ntermd_1(i,j,k,iblock)
1204 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1205 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1206 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1207 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1210 write (iout,*) 'Pairs of angles:'
1211 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1212 do l=1,ntermd_2(i,j,k,iblock)
1213 write (iout,'(i5,20f10.5)') &
1214 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1217 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1218 do l=1,ntermd_2(i,j,k,iblock)
1219 write (iout,'(i5,20f10.5)') &
1220 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1221 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1230 !elwrite(iout,*) "parmread kontrol sc-bb"
1231 ! Read of Side-chain backbone correlation parameters
1232 ! Modified 11 May 2012 by Adasko
1235 read (isccor,*) nsccortyp
1238 !c maxinter is maximum interaction sites
1239 !write(iout,*)"maxterm_sccor",maxterm_sccor
1240 !el from module energy-------------
1241 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1242 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1243 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1244 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1245 !-----------------------------------
1246 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1247 !-----------------------------------
1248 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1249 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1250 -nsccortyp:nsccortyp))
1251 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1252 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1253 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1254 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1255 !-----------------------------------
1256 do i=-nsccortyp,nsccortyp
1257 do j=-nsccortyp,nsccortyp
1261 !-----------------------------------
1263 read (isccor,*) (isccortyp(i),i=1,ntyp)
1265 isccortyp(i)=-isccortyp(-i)
1267 iscprol=isccortyp(20)
1268 ! write (iout,*) 'ntortyp',ntortyp
1270 !c maxinter is maximum interaction sites
1275 nterm_sccor(i,j),nlor_sccor(i,j)
1281 nterm_sccor(-i,j)=nterm_sccor(i,j)
1282 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1283 nterm_sccor(i,-j)=nterm_sccor(i,j)
1284 do k=1,nterm_sccor(i,j)
1285 read (isccor,*) kk,v1sccor(k,l,i,j),&
1287 if (j.eq.iscprol) then
1288 if (i.eq.isccortyp(10)) then
1289 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1290 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1292 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1293 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1294 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1295 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1296 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1297 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1298 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1299 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1302 if (i.eq.isccortyp(10)) then
1303 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1304 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1306 if (j.eq.isccortyp(10)) then
1307 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1308 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1310 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1311 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1312 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1313 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1314 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1315 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1319 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1320 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1321 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1322 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1325 do k=1,nlor_sccor(i,j)
1326 read (isccor,*) kk,vlor1sccor(k,i,j),&
1327 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1328 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1329 (1+vlor3sccor(k,i,j)**2)
1331 v0sccor(l,i,j)=v0ijsccor
1332 v0sccor(l,-i,j)=v0ijsccor1
1333 v0sccor(l,i,-j)=v0ijsccor2
1334 v0sccor(l,-i,-j)=v0ijsccor3
1340 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1343 write (iout,*) 'ityp',i,' jtyp',j
1344 write (iout,*) 'Fourier constants'
1345 do k=1,nterm_sccor(i,j)
1346 write (iout,'(2(1pe15.5))') &
1347 (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
1349 write (iout,*) 'Lorenz constants'
1350 do k=1,nlor_sccor(i,j)
1351 write (iout,'(3(1pe15.5))') &
1352 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1358 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1359 ! interaction energy of the Gly, Ala, and Pro prototypes.
1361 read (ifourier,*) nloctyp
1362 !el write(iout,*)"nloctyp",nloctyp
1363 !el from module energy-------
1364 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1365 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1366 allocate(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1367 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1368 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1369 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1370 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1371 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1373 do ii=-nloctyp-1,nloctyp+1
1381 ctilde(j,i,ii)=0.0d0
1382 dtilde(j,i,ii)=0.0d0
1386 !--------------------------------
1387 allocate(b(13,0:nloctyp))
1391 read (ifourier,*) (b(ii,i),ii=1,13)
1393 write (iout,*) 'Type',i
1394 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1402 B1tilde(1,i) = b(3,i)
1403 B1tilde(2,i) =-b(5,i)
1404 B1tilde(1,-i) =-b(3,i)
1405 B1tilde(2,-i) =b(5,i)
1406 ! b1tilde(1,i)=0.0d0
1407 ! b1tilde(2,i)=0.0d0
1427 Ctilde(1,1,i)=b(7,i)
1428 Ctilde(1,2,i)=b(9,i)
1429 Ctilde(2,1,i)=-b(9,i)
1430 Ctilde(2,2,i)=b(7,i)
1431 Ctilde(1,1,-i)=b(7,i)
1432 Ctilde(1,2,-i)=-b(9,i)
1433 Ctilde(2,1,-i)=b(9,i)
1434 Ctilde(2,2,-i)=b(7,i)
1436 ! Ctilde(1,1,i)=0.0d0
1437 ! Ctilde(1,2,i)=0.0d0
1438 ! Ctilde(2,1,i)=0.0d0
1439 ! Ctilde(2,2,i)=0.0d0
1452 Dtilde(1,1,i)=b(6,i)
1453 Dtilde(1,2,i)=b(8,i)
1454 Dtilde(2,1,i)=-b(8,i)
1455 Dtilde(2,2,i)=b(6,i)
1456 Dtilde(1,1,-i)=b(6,i)
1457 Dtilde(1,2,-i)=-b(8,i)
1458 Dtilde(2,1,-i)=b(8,i)
1459 Dtilde(2,2,-i)=b(6,i)
1461 ! Dtilde(1,1,i)=0.0d0
1462 ! Dtilde(1,2,i)=0.0d0
1463 ! Dtilde(2,1,i)=0.0d0
1464 ! Dtilde(2,2,i)=0.0d0
1465 EE(1,1,i)= b(10,i)+b(11,i)
1466 EE(2,2,i)=-b(10,i)+b(11,i)
1467 EE(2,1,i)= b(12,i)-b(13,i)
1468 EE(1,2,i)= b(12,i)+b(13,i)
1469 EE(1,1,-i)= b(10,i)+b(11,i)
1470 EE(2,2,-i)=-b(10,i)+b(11,i)
1471 EE(2,1,-i)=-b(12,i)+b(13,i)
1472 EE(1,2,-i)=-b(12,i)-b(13,i)
1478 ! ee(2,1,i)=ee(1,2,i)
1483 write (iout,*) 'Type',i
1485 ! write (iout,'(f10.5)') B1(:,i)
1486 write(iout,*) B1(1,i),B1(2,i)
1488 ! write (iout,'(f10.5)') B2(:,i)
1489 write(iout,*) B2(1,i),B2(2,i)
1492 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1496 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1500 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1505 ! Read electrostatic-interaction parameters
1508 write (iout,'(/a)') 'Electrostatic interaction constants:'
1509 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1510 'IT','JT','APP','BPP','AEL6','AEL3'
1512 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1513 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1514 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1515 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1520 app (i,j)=epp(i,j)*rri*rri
1521 bpp (i,j)=-2.0D0*epp(i,j)*rri
1522 ael6(i,j)=elpp6(i,j)*4.2D0**6
1523 ael3(i,j)=elpp3(i,j)*4.2D0**3
1524 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1529 ! Read side-chain interaction parameters.
1531 !el from module energy - COMMON.INTERACT-------
1532 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1533 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1534 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1535 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1536 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1547 !--------------------------------
1549 read (isidep,*) ipot,expon
1550 !el if (ipot.lt.1 .or. ipot.gt.5) then
1551 ! write (iout,'(2a)') 'Error while reading SC interaction',&
1552 ! 'potential file - unknown potential type.'
1556 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1557 ', exponents are ',expon,2*expon
1558 ! goto (10,20,30,30,40) ipot
1560 !----------------------- LJ potential ---------------------------------
1562 ! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1563 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1565 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1566 write (iout,'(a/)') 'The epsilon array:'
1567 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1568 write (iout,'(/a)') 'One-body parameters:'
1569 write (iout,'(a,4x,a)') 'residue','sigma'
1570 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1573 !----------------------- LJK potential --------------------------------
1575 ! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1576 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1577 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1579 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1580 write (iout,'(a/)') 'The epsilon array:'
1581 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1582 write (iout,'(/a)') 'One-body parameters:'
1583 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1584 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1588 !---------------------- GB or BP potential -----------------------------
1592 read (isidep,*)(eps(i,j),j=i,ntyp)
1594 read (isidep,*)(sigma0(i),i=1,ntyp)
1595 read (isidep,*)(sigii(i),i=1,ntyp)
1596 read (isidep,*)(chip(i),i=1,ntyp)
1597 read (isidep,*)(alp(i),i=1,ntyp)
1598 ! For the GB potential convert sigma'**2 into chi'
1601 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1605 write (iout,'(/a/)') 'Parameters of the BP potential:'
1606 write (iout,'(a/)') 'The epsilon array:'
1607 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1608 write (iout,'(/a)') 'One-body parameters:'
1609 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
1611 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
1612 chip(i),alp(i),i=1,ntyp)
1615 !--------------------- GBV potential -----------------------------------
1617 ! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1618 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1619 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
1620 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1622 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1623 write (iout,'(a/)') 'The epsilon array:'
1624 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1625 write (iout,'(/a)') 'One-body parameters:'
1626 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
1627 's||/s_|_^2',' chip ',' alph '
1628 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
1629 sigii(i),chip(i),alp(i),i=1,ntyp)
1632 write (iout,'(2a)') 'Error while reading SC interaction',&
1633 'potential file - unknown potential type.'
1639 !-----------------------------------------------------------------------
1640 ! Calculate the "working" parameters of SC interactions.
1642 !el from module energy - COMMON.INTERACT-------
1643 allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
1644 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
1654 !--------------------------------
1663 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1664 sigma(j,i)=sigma(i,j)
1665 rs0(i,j)=dwa16*sigma(i,j)
1669 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
1670 'Working parameters of the SC interactions:',&
1671 ' a ',' b ',' augm ',' sigma ',' r0 ',&
1676 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1685 sigeps=dsign(1.0D0,epsij)
1687 aa(i,j)=epsij*rrij*rrij
1688 bb(i,j)=-sigeps*epsij*rrij
1692 sigt1sq=sigma0(i)**2
1693 sigt2sq=sigma0(j)**2
1696 ratsig1=sigt2sq/sigt1sq
1697 ratsig2=1.0D0/ratsig1
1698 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1699 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1700 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1704 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1705 sigmaii(i,j)=rsum_max
1706 sigmaii(j,i)=rsum_max
1708 ! sigmaii(i,j)=r0(i,j)
1709 ! sigmaii(j,i)=r0(i,j)
1711 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1712 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1713 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1714 augm(i,j)=epsij*r_augm**(2*expon)
1715 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1722 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
1723 restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
1724 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1729 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
1737 ! Define the SC-p interaction constants
1747 !elwrite(iout,*) "parmread kontrol before oldscp"
1749 ! Define the SC-p interaction constants
1753 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1755 ! aad(i,1)=0.3D0*4.0D0**12
1756 ! Following line for constants currently implemented
1757 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
1758 aad(i,1)=1.5D0*4.0D0**12
1759 ! aad(i,1)=0.17D0*5.6D0**12
1761 ! "Soft" SC-p repulsion
1763 ! Following line for constants currently implemented
1764 ! aad(i,1)=0.3D0*4.0D0**6
1765 ! "Hard" SC-p repulsion
1766 bad(i,1)=3.0D0*4.0D0**6
1767 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1776 ! 8/9/01 Read the SC-p interaction constants from file
1779 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1782 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1783 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1784 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1785 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1789 write (iout,*) "Parameters of SC-p interactions:"
1791 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
1792 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1797 ! Define the constants of the disulfide bridge
1801 ! Old arbitrary potential - commented out.
1806 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1807 ! energy surface of diethyl disulfide.
1808 ! A. Liwo and U. Kozlowska, 11/24/03
1819 write (iout,'(/a)') "Disulfide bridge parameters:"
1820 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1821 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1822 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1823 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
1827 end subroutine parmread
1829 !-----------------------------------------------------------------------------
1831 !-----------------------------------------------------------------------------
1832 subroutine mygetenv(string,var)
1836 ! This subroutine passes the environmental variables to FORTRAN program.
1837 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
1838 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
1839 ! reads the environmental variables from $HOME/.env
1841 ! Usage: As for the standard FORTRAN GETENV subroutine.
1843 ! Purpose: some versions/installations of MPI do not transfer the environmental
1844 ! variables to slave processors, if these variables are set in the shell script
1845 ! from which mpirun is called.
1854 character*(*) :: string,var
1855 #if defined(MYGETENV) && defined(MPI)
1856 ! include "DIMENSIONS.ZSCOPT"
1858 ! include "COMMON.MPI"
1859 !el character*360 ucase
1861 character(len=360) :: string1(360),karta
1862 character(len=240) :: home
1865 call getenv("HOME",home)
1866 open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
1868 read (99,end=111,err=111,'(a)') karta
1872 call split_string(karta,string1,80,n)
1873 if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
1874 string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
1876 print *,"Processor",me,": ",var(:ilen(var)),&
1877 " assigned to ",string(:ilen(string))
1882 111 print *,"Environment variable ",string(:ilen(string))," not set."
1885 112 print *,"Error opening environment file!"
1887 call getenv(string,var)
1890 end subroutine mygetenv
1891 !-----------------------------------------------------------------------------
1893 !-----------------------------------------------------------------------------
1894 subroutine read_general_data(*)
1896 use control_data, only:indpdb,symetr
1897 use energy_data, only:distchainmax
1899 ! include "DIMENSIONS"
1900 ! include "DIMENSIONS.ZSCOPT"
1901 ! include "DIMENSIONS.FREE"
1902 ! include "COMMON.TORSION"
1903 ! include "COMMON.INTERACT"
1904 ! include "COMMON.IOUNITS"
1905 ! include "COMMON.TIME1"
1906 ! include "COMMON.PROT"
1907 ! include "COMMON.PROTFILES"
1908 ! include "COMMON.CHAIN"
1909 ! include "COMMON.NAMES"
1910 ! include "COMMON.FFIELD"
1911 ! include "COMMON.ENEPS"
1912 ! include "COMMON.WEIGHTS"
1913 ! include "COMMON.FREE"
1914 ! include "COMMON.CONTROL"
1915 ! include "COMMON.ENERGIES"
1916 character(len=800) :: controlcard
1917 integer :: i,j,k,ii,n_ene_found
1918 integer :: ind,itype1,itype2,itypf,itypsc,itypp
1921 !el character*16 ucase
1922 character(len=16) :: key
1924 call card_concat(controlcard,.true.)
1925 call readi(controlcard,"N_ENE",n_eneW,max_eneW)
1926 if (n_eneW.gt.max_eneW) then
1927 write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
1931 call readi(controlcard,"NPARMSET",nparmset,1)
1932 !elwrite(iout,*)"in read_gen data"
1933 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
1934 call readi(controlcard,"IPARMPRINT",iparmprint,1)
1935 write (iout,*) "PARMPRINT",iparmprint
1936 if (nparmset.gt.max_parm) then
1937 write (iout,*) "Error: parameter out of range: NPARMSET",&
1941 !elwrite(iout,*)"in read_gen data"
1942 call readi(controlcard,"MAXIT",maxit,5000)
1943 call reada(controlcard,"FIMIN",fimin,1.0d-3)
1944 call readi(controlcard,"ENSEMBLES",ensembles,0)
1945 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
1946 write (iout,*) "Number of energy parameter sets",nparmset
1947 allocate(isampl(nparmset))
1948 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
1949 write (iout,*) "MaxSlice",MaxSlice
1950 call readi(controlcard,"NSLICE",nslice,1)
1951 !elwrite(iout,*)"in read_gen data"
1953 if (nslice.gt.MaxSlice) then
1954 write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
1958 write (iout,*) "Frequency of storing conformations",&
1959 (isampl(i),i=1,nparmset)
1960 write (iout,*) "Maxit",maxit," Fimin",fimin
1961 call readi(controlcard,"NQ",nQ,1)
1962 if (nQ.gt.MaxQ) then
1963 write (iout,*) "Error: parameter out of range: NQ",nq,&
1968 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
1969 call reada(controlcard,"DELTA",delta,1.0d-2)
1970 call readi(controlcard,"EINICHECK",einicheck,2)
1971 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
1972 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
1973 call readi(controlcard,"RESCALE",rescale_modeW,1)
1974 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
1975 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1976 call readi(controlcard,'SYM',symetr,1)
1977 write (iout,*) "DISTCHAINMAX",distchainmax
1978 write (iout,*) "delta",delta
1979 write (iout,*) "einicheck",einicheck
1980 write (iout,*) "rescale_mode",rescale_modeW
1982 bxfile=index(controlcard,"BXFILE").gt.0
1983 cxfile=index(controlcard,"CXFILE").gt.0
1984 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
1986 histfile=index(controlcard,"HISTFILE").gt.0
1987 histout=index(controlcard,"HISTOUT").gt.0
1988 entfile=index(controlcard,"ENTFILE").gt.0
1989 zscfile=index(controlcard,"ZSCFILE").gt.0
1990 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1991 write (iout,*) "with_dihed_constr ",with_dihed_constr
1992 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
1994 end subroutine read_general_data
1995 !------------------------------------------------------------------------------
1996 subroutine read_efree(*)
1998 ! Read molecular data
2001 ! include 'DIMENSIONS'
2002 ! include 'DIMENSIONS.ZSCOPT'
2003 ! include 'DIMENSIONS.COMPAR'
2004 ! include 'DIMENSIONS.FREE'
2005 ! include 'COMMON.IOUNITS'
2006 ! include 'COMMON.TIME1'
2007 ! include 'COMMON.SBRIDGE'
2008 ! include 'COMMON.CONTROL'
2009 ! include 'COMMON.CHAIN'
2010 ! include 'COMMON.HEADER'
2011 ! include 'COMMON.GEO'
2012 ! include 'COMMON.FREE'
2013 character(len=320) :: controlcard !,ucase
2014 integer :: iparm,ib,i,j,npars
2024 ! call alloc_wham_arrays
2025 ! allocate(nT_h(nParmSet))
2026 ! allocate(replica(nParmSet))
2027 ! allocate(umbrella(nParmSet))
2028 ! allocate(read_iset(nParmSet))
2029 ! allocate(nT_h(nParmSet))
2033 call card_concat(controlcard,.true.)
2034 call readi(controlcard,'NT',nT_h(iparm),1)
2035 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
2037 if (nT_h(iparm).gt.MaxT_h) then
2038 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
2042 replica(iparm)=index(controlcard,"REPLICA").gt.0
2043 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
2044 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
2045 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
2046 replica(iparm)," umbrella ",umbrella(iparm),&
2047 " read_iset",read_iset(iparm)
2050 call card_concat(controlcard,.true.)
2051 call readi(controlcard,'NR',nR(ib,iparm),1)
2052 if (umbrella(iparm)) then
2055 nRR(ib,iparm)=nR(ib,iparm)
2057 if (nR(ib,iparm).gt.MaxR) then
2058 write (iout,*) "Error: parameter out of range: NR",&
2062 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
2063 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
2064 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
2067 call card_concat(controlcard,.true.)
2068 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
2070 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
2075 write (iout,*) "ib",ib," beta_h",&
2076 1.0d0/(0.001987*beta_h(ib,iparm))
2077 write (iout,*) "nR",nR(ib,iparm)
2078 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
2080 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
2081 "q0",(q0(j,i,ib,iparm),j=1,nQ)
2093 nR(ib,iparm)=nR(ib,1)
2094 if (umbrella(iparm)) then
2097 nRR(ib,iparm)=nR(ib,1)
2099 beta_h(ib,iparm)=beta_h(ib,1)
2101 f(i,ib,iparm)=f(i,ib,1)
2103 KH(j,i,ib,iparm)=KH(j,i,ib,1)
2104 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
2107 replica(iparm)=replica(1)
2108 umbrella(iparm)=umbrella(1)
2109 read_iset(iparm)=read_iset(1)
2116 end subroutine read_efree
2117 !-----------------------------------------------------------------------------
2118 subroutine read_protein_data(*)
2120 ! include "DIMENSIONS"
2121 ! include "DIMENSIONS.ZSCOPT"
2122 ! include "DIMENSIONS.FREE"
2126 integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
2127 ! include "COMMON.MPI"
2129 ! include "COMMON.CHAIN"
2130 ! include "COMMON.IOUNITS"
2131 ! include "COMMON.PROT"
2132 ! include "COMMON.PROTFILES"
2133 ! include "COMMON.NAMES"
2134 ! include "COMMON.FREE"
2135 ! include "COMMON.OBCINKA"
2136 character(len=64) :: nazwa
2137 character(len=16000) :: controlcard
2138 integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
2139 !el external ilen,iroof
2148 ! Read names of files with conformation data.
2149 if (replica(iparm)) then
2155 do ii=1,nRR(ib,iparm)
2156 write (iout,*) "Parameter set",iparm," temperature",ib,&
2159 call card_concat(controlcard,.true.)
2160 write (iout,*) controlcard(:ilen(controlcard))
2161 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
2162 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
2163 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
2164 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
2165 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
2166 maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
2167 call reada(controlcard,"TIME_START",&
2168 time_start_collect(ii,ib,iparm),0.0d0)
2169 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
2171 write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
2172 " rec_end",rec_end(ii,ib,iparm)
2173 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
2174 " time_end",time_end_collect(ii,ib,iparm)
2176 if (replica(iparm)) then
2177 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
2178 write (iout,*) "Number of trajectories",totraj(ii,iparm)
2181 if (nfile_bin(ii,ib,iparm).lt.2 &
2182 .and. nfile_asc(ii,ib,iparm).eq.0 &
2183 .and. nfile_cx(ii,ib,iparm).eq.0) then
2184 write (iout,*) "Error - no action specified!"
2187 if (nfile_bin(ii,ib,iparm).gt.0) then
2188 call card_concat(controlcard,.false.)
2189 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
2190 maxfile_prot,nfile_bin(ii,ib,iparm))
2192 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
2193 write(iout,*) (protfiles(i,1,ii,ib,iparm),&
2194 i=1,nfile_bin(ii,ib,iparm))
2197 if (nfile_asc(ii,ib,iparm).gt.0) then
2198 call card_concat(controlcard,.false.)
2199 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2200 maxfile_prot,nfile_asc(ii,ib,iparm))
2202 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
2203 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2204 i=1,nfile_asc(ii,ib,iparm))
2206 else if (nfile_cx(ii,ib,iparm).gt.0) then
2207 call card_concat(controlcard,.false.)
2208 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2209 maxfile_prot,nfile_cx(ii,ib,iparm))
2211 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
2212 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2213 i=1,nfile_cx(ii,ib,iparm))
2222 end subroutine read_protein_data
2223 !-------------------------------------------------------------------------------
2224 subroutine readsss(rekord,lancuch,wartosc,default)
2226 character*(*) :: rekord,lancuch,wartosc,default
2227 character(len=80) :: aux
2228 integer :: lenlan,lenrec,iread,ireade
2232 lenlan=ilen(lancuch)
2234 iread=index(rekord,lancuch(:lenlan)//"=")
2235 ! print *,"rekord",rekord," lancuch",lancuch
2236 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2237 if (iread.eq.0) then
2241 iread=iread+lenlan+1
2242 ! print *,"iread",iread
2243 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2244 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2246 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2248 ! print *,"iread",iread
2249 if (iread.gt.lenrec) then
2254 ! print *,"ireade",ireade
2255 do while (ireade.lt.lenrec .and. &
2256 .not.iblnk(rekord(ireade:ireade)))
2259 wartosc=rekord(iread:ireade)
2261 end subroutine readsss
2262 !----------------------------------------------------------------------------
2263 subroutine multreads(rekord,lancuch,tablica,dim,default)
2266 character*(*) rekord,lancuch,tablica(dim),default
2267 character(len=80) :: aux
2268 integer :: lenlan,lenrec,iread,ireade
2275 lenlan=ilen(lancuch)
2277 iread=index(rekord,lancuch(:lenlan)//"=")
2278 ! print *,"rekord",rekord," lancuch",lancuch
2279 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2280 if (iread.eq.0) return
2281 iread=iread+lenlan+1
2283 ! print *,"iread",iread
2284 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2285 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2287 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2289 ! print *,"iread",iread
2290 if (iread.gt.lenrec) return
2292 ! print *,"ireade",ireade
2293 do while (ireade.lt.lenrec .and. &
2294 .not.iblnk(rekord(ireade:ireade)))
2297 tablica(i)=rekord(iread:ireade)
2300 end subroutine multreads
2301 !----------------------------------------------------------------------------
2302 subroutine split_string(rekord,tablica,dim,nsub)
2304 integer :: dim,nsub,i,ii,ll,kk
2305 character*(*) tablica(dim)
2306 character*(*) rekord
2316 ! Find the start of term name
2318 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2321 ! Parse the name into TABLICA(i) until blank found
2322 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2324 tablica(i)(kk:kk)=rekord(ii:ii)
2327 if (kk.gt.0) nsub=nsub+1
2328 if (ii.gt.ll) return
2331 end subroutine split_string
2332 !--------------------------------------------------------------------------------
2334 !--------------------------------------------------------------------------------
2335 subroutine read_compar
2337 ! Read molecular data
2339 use conform_compar, only:alloc_compar_arrays
2340 use control_data, only:pdbref
2341 use geometry_data, only:deg2rad,rad2deg
2343 ! include 'DIMENSIONS'
2344 ! include 'DIMENSIONS.ZSCOPT'
2345 ! include 'DIMENSIONS.COMPAR'
2346 ! include 'DIMENSIONS.FREE'
2347 ! include 'COMMON.IOUNITS'
2348 ! include 'COMMON.TIME1'
2349 ! include 'COMMON.SBRIDGE'
2350 ! include 'COMMON.CONTROL'
2351 ! include 'COMMON.COMPAR'
2352 ! include 'COMMON.CHAIN'
2353 ! include 'COMMON.HEADER'
2354 ! include 'COMMON.GEO'
2355 ! include 'COMMON.FREE'
2356 character(len=320) :: controlcard !,ucase
2357 character(len=64) :: wfile
2361 !elwrite(iout,*)"jestesmy w read_compar"
2362 call card_concat(controlcard,.true.)
2363 pdbref=(index(controlcard,'PDBREF').gt.0)
2364 call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
2365 call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
2366 call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
2367 call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
2368 verbose = index(controlcard,"VERBOSE").gt.0
2369 lgrp=index(controlcard,"STATIN").gt.0
2370 lgrp_out=index(controlcard,"STATOUT").gt.0
2371 merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
2372 binary = index(controlcard,"BINARY").gt.0
2373 rmscut_base_up=rmscut_base_up/50
2374 rmscut_base_low=rmscut_base_low/50
2375 call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
2376 call readi(controlcard,'NLEVEL',nlevel,1)
2377 if (nlevel.lt.0) then
2379 call alloc_compar_arrays(maxfrag,1)
2382 allocate(nfrag(nlevel))
2384 ! Read the data pertaining to elementary fragments (level 1)
2385 call readi(controlcard,'NFRAG',nfrag(1),0)
2386 write(iout,*)"nfrag(1)",nfrag(1)
2387 call alloc_compar_arrays(nfrag(1),nlevel)
2389 call card_concat(controlcard,.true.)
2390 write (iout,*) controlcard(:ilen(controlcard))
2391 call readi(controlcard,'NPIECE',npiece(j,1),0)
2392 call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
2393 call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
2394 call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
2395 call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
2396 call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
2397 call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
2398 call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
2399 call readi(controlcard,'RMS',irms(j,1),0)
2400 call readi(controlcard,'LOCAL',iloc(j),1)
2401 call readi(controlcard,'ELCONT',ielecont(j,1),1)
2402 if (ielecont(j,1).eq.0) then
2403 call readi(controlcard,'SCCONT',isccont(j,1),1)
2405 ang_cut(j)=ang_cut(j)*deg2rad
2406 ang_cut1(j)=ang_cut1(j)*deg2rad
2408 call card_concat(controlcard,.true.)
2409 call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
2410 call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
2412 write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
2413 (ifrag(1,k,j),ifrag(2,k,j),&
2414 k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
2415 " ang_cut1",ang_cut1(j)*rad2deg
2416 write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
2417 write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
2418 write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
2419 " ilocal",iloc(j)," isccont",isccont(j,1)
2421 ! Read data pertaning to higher levels
2423 call card_concat(controlcard,.true.)
2424 call readi(controlcard,'NFRAG',NFRAG(i),0)
2425 write (iout,*) "i",i," nfrag",nfrag(i)
2427 call card_concat(controlcard,.true.)
2429 call readi(controlcard,'ELCONT',ielecont(j,i),0)
2430 if (ielecont(j,i).eq.0) then
2431 call readi(controlcard,'SCCONT',isccont(j,i),1)
2433 call readi(controlcard,'RMS',irms(j,i),0)
2439 call readi(controlcard,'NPIECE',npiece(j,i),0)
2440 call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
2441 call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
2442 call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
2444 call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
2445 call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
2446 write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
2447 n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
2448 " isccont",isccont(j,i)," irms",irms(j,i)
2449 write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
2450 write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
2451 write(iout,*)"nc_frac",nc_fragm(j,i),&
2452 " nc_req",nc_req_setf(j,i)
2455 if (binary) write (iout,*) "Classes written in binary format."
2458 call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
2459 call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
2460 call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
2461 call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
2462 call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
2463 call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
2464 call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
2465 call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
2466 call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
2467 call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
2468 call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
2469 call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
2470 call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
2471 call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
2472 call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
2473 call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
2474 call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
2475 call readi(controlcard,'RMS_SINGLE',irms_single,0)
2476 call readi(controlcard,'CONT_SINGLE',icont_single,1)
2477 call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
2478 call readi(controlcard,'RMS_PAIR',irms_pair,0)
2479 call readi(controlcard,'CONT_PAIR',icont_pair,1)
2480 call readi(controlcard,'SPLIT_BET',isplit_bet,0)
2481 angcut_hel=angcut_hel*deg2rad
2482 angcut1_hel=angcut1_hel*deg2rad
2483 angcut_bet=angcut_bet*deg2rad
2484 angcut1_bet=angcut1_bet*deg2rad
2485 angcut_strand=angcut_strand*deg2rad
2486 angcut1_strand=angcut1_strand*deg2rad
2487 write (iout,*) "Automatic detection of structural elements"
2488 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2489 ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2490 ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
2491 ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2492 ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
2493 ' SPLIT_BET',isplit_bet
2494 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2495 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2496 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2497 ' MAXANG_HEL',angcut1_hel*rad2deg
2498 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2499 ' MAXANG_BET',angcut1_bet*rad2deg
2500 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2501 ' MAXANG_STRAND',angcut1_strand*rad2deg
2502 write (iout,*) 'FRAC_MIN',frac_min_set
2504 end subroutine read_compar
2505 !--------------------------------------------------------------------------------
2507 !--------------------------------------------------------------------------------
2508 subroutine read_ref_structure(*)
2510 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
2513 use control_data, only:pdbref
2514 use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
2515 nstart_sup,nstart_seq,nperm,nres0
2516 use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype
2517 use compare, only:seq_comp !,contact,elecont
2518 use geometry, only:chainbuild,dist
2519 use io_config, only:readpdb
2521 use conform_compar, only:contact,elecont
2523 ! include 'DIMENSIONS'
2524 ! include 'DIMENSIONS.ZSCOPT'
2525 ! include 'DIMENSIONS.COMPAR'
2526 ! include 'COMMON.IOUNITS'
2527 ! include 'COMMON.GEO'
2528 ! include 'COMMON.VAR'
2529 ! include 'COMMON.INTERACT'
2530 ! include 'COMMON.LOCAL'
2531 ! include 'COMMON.NAMES'
2532 ! include 'COMMON.CHAIN'
2533 ! include 'COMMON.FFIELD'
2534 ! include 'COMMON.SBRIDGE'
2535 ! include 'COMMON.HEADER'
2536 ! include 'COMMON.CONTROL'
2537 ! include 'COMMON.CONTACTS1'
2538 ! include 'COMMON.PEPTCONT'
2539 ! include 'COMMON.TIME1'
2540 ! include 'COMMON.COMPAR'
2541 character(len=4) :: sequence(nres)
2543 !el real(kind=8) :: x(maxvar)
2544 integer :: itype_pdb(nres)
2545 !el logical seq_comp
2546 integer :: i,j,k,nres_pdb,iaux
2547 real(kind=8) :: ddsc !el,dist
2548 integer :: kkk !,ilen
2552 write (iout,*) "pdbref",pdbref
2554 read(inp,'(a)') pdbfile
2555 write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
2556 pdbfile(:ilen(pdbfile))
2557 open(ipdbin,file=pdbfile,status='old',err=33)
2559 33 write (iout,'(a)') 'Error opening PDB file.'
2563 itype_pdb(i)=itype(i)
2570 itype_pdb(i)=itype(i)
2578 if (nsup.le.(nct-nnt+1)) then
2579 do i=0,nct-nnt+1-nsup
2580 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),&
2582 do j=nnt+nsup-1,nnt,-1
2584 cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
2587 do j=nnt+nsup-1,nnt,-1
2589 cref(k,j+i,kkk)=cref(k,j,kkk)
2591 phi_ref(j+i)=phi_ref(j)
2592 theta_ref(j+i)=theta_ref(j)
2593 alph_ref(j+i)=alph_ref(j)
2594 omeg_ref(j+i)=omeg_ref(j)
2598 write (iout,'(i5,3f10.5,5x,3f10.5)') &
2599 j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
2607 write (iout,'(a)') &
2608 'Error - sequences to be superposed do not match.'
2611 do i=0,nsup-(nct-nnt+1)
2612 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),&
2615 nstart_sup=nstart_sup+i
2620 write (iout,'(a)') &
2621 'Error - sequences to be superposed do not match.'
2625 write (iout,'(a,i5)') &
2626 'Experimental structure begins at residue',nstart_seq
2628 call read_angles(inp,*38)
2630 38 write (iout,'(a)') 'Error reading reference structure.'
2639 cref(j,i,kkk)=c(j,i)
2643 nend_sup=nstart_sup+nsup-1
2646 c(j,i)=cref(j,i,kkk)
2651 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
2653 if (itype(i).ne.10) then
2654 ddsc = dist(i,nres+i)
2656 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
2660 dc_norm(j,nres+i)=0.0d0
2663 ! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
2664 ! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
2665 ! dc_norm(3,nres+i)**2
2667 dc(j,i)=c(j,i+1)-c(j,i)
2671 dc_norm(j,i)=dc(j,i)/ddsc
2674 ! print *,"Calling contact"
2675 call contact(.true.,ncont_ref,icont_ref(1,1),&
2676 nstart_sup,nend_sup)
2677 ! print *,"Calling elecont"
2678 call elecont(.true.,ncont_pept_ref,&
2679 icont_pept_ref(1,1),&
2680 nstart_sup,nend_sup)
2681 write (iout,'(a,i3,a,i3,a,i3,a)') &
2682 'Number of residues to be superposed:',nsup,&
2683 ' (from residue',nstart_sup,' to residue',&
2686 end subroutine read_ref_structure
2687 !--------------------------------------------------------------------------------
2689 !--------------------------------------------------------------------------------
2690 subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
2692 use geometry_data, only:nres,c
2693 use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype
2694 ! implicit real*8 (a-h,o-z)
2695 ! include 'DIMENSIONS'
2696 ! include 'DIMENSIONS.ZSCOPT'
2697 ! include 'COMMON.CHAIN'
2698 ! include 'COMMON.INTERACT'
2699 ! include 'COMMON.NAMES'
2700 ! include 'COMMON.IOUNITS'
2701 ! include 'COMMON.HEADER'
2702 ! include 'COMMON.SBRIDGE'
2703 character(len=50) :: tytul
2704 character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
2705 'D','E','F','G','H','I','J'/),shape(chainid))
2706 integer,dimension(nres) :: ica !(maxres)
2707 real(kind=8) :: temp,efree,etot,entropy,rmsdev
2708 integer :: ii,i,j,iti,ires,iatom,ichain
2709 write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
2711 write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
2713 write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
2720 if (iti.eq.ntyp1) then
2723 write (ipdb,'(a)') 'TER'
2728 write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
2732 write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
2733 ires,(c(j,nres+i),j=1,3)
2737 write (ipdb,'(a)') 'TER'
2739 if (itype(i).eq.ntyp1) cycle
2740 if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
2741 write (ipdb,30) ica(i),ica(i+1)
2742 else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
2743 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2744 else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
2745 write (ipdb,30) ica(i),ica(i)+1
2748 if (itype(nct).ne.10) then
2749 write (ipdb,30) ica(nct),ica(nct)+1
2752 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2754 write (ipdb,'(a6)') 'ENDMDL'
2755 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2756 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2757 30 FORMAT ('CONECT',8I5)
2759 end subroutine pdboutW
2761 !------------------------------------------------------------------------------
2763 !-----------------------------------------------------------------------------
2764 !-----------------------------------------------------------------------------