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: maxterm,maxlor,maxterm_sccor,& !maxtor
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(-nthetyp-1:nthetyp+1,&
687 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
688 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
689 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
690 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
691 !(maxtheterm,-maxthetyp1:maxthetyp1,&
692 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
693 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
694 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
695 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
696 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
697 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
698 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
699 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
700 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
701 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
702 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
703 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
704 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
705 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
706 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,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"
716 aa0thet(:,:,:,:)=0.0d0
717 aathet(:,:,:,:,:)=0.0d0
718 bbthet(:,:,:,:,:,:)=0.0d0
719 ccthet(:,:,:,:,:,:)=0.0d0
720 ddthet(:,:,:,:,:,:)=0.0d0
721 eethet(:,:,:,:,:,:)=0.0d0
722 ffthet(:,:,:,:,:,:,:)=0.0d0
723 ggthet(:,:,:,:,:,:,:)=0.0d0
727 do j=-nthetyp,nthetyp
728 do k=-nthetyp,nthetyp
729 read (ithep,'(6a)') res1
730 read (ithep,*) aa0thet(i,j,k,iblock)
731 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
733 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
734 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
735 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
736 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
739 (((ffthet(llll,lll,ll,i,j,k,iblock),&
740 ffthet(lll,llll,ll,i,j,k,iblock),&
741 ggthet(llll,lll,ll,i,j,k,iblock),&
742 ggthet(lll,llll,ll,i,j,k,iblock),&
743 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
748 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
749 ! coefficients of theta-and-gamma-dependent terms are zero.
754 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
755 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
757 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
758 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
761 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
763 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
766 ! Substitution for D aminoacids from symmetry.
769 do j=-nthetyp,nthetyp
770 do k=-nthetyp,nthetyp
771 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
773 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
777 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
778 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
779 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
780 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
786 ffthet(llll,lll,ll,i,j,k,iblock)= &
787 ffthet(llll,lll,ll,-i,-j,-k,iblock)
788 ffthet(lll,llll,ll,i,j,k,iblock)= &
789 ffthet(lll,llll,ll,-i,-j,-k,iblock)
790 ggthet(llll,lll,ll,i,j,k,iblock)= &
791 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
792 ggthet(lll,llll,ll,i,j,k,iblock)= &
793 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
803 ! Control printout of the coefficients of virtual-bond-angle potentials
807 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
811 write (iout,'(//4a)') &
812 'Type ',onelett(i),onelett(j),onelett(k)
813 write (iout,'(//a,10x,a)') " l","a[l]"
814 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
815 write (iout,'(i2,1pe15.5)') &
816 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
818 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
819 "b",l,"c",l,"d",l,"e",l
821 write (iout,'(i2,4(1pe15.5))') m,&
822 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
823 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
827 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
828 "f+",l,"f-",l,"g+",l,"g-",l
831 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
832 ffthet(n,m,l,i,j,k,iblock),&
833 ffthet(m,n,l,i,j,k,iblock),&
834 ggthet(n,m,l,i,j,k,iblock),&
835 ggthet(m,n,l,i,j,k,iblock)
846 !-------------------------------------------
847 allocate(nlob(ntyp1)) !(ntyp1)
848 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
849 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
850 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
868 gaussc(l,k,j,i)=0.0D0
876 ! Read the parameters of the probability distribution/energy expression
877 ! of the side chains.
880 !c write (iout,*) "tu dochodze",i
881 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
885 dsc_inv(i)=1.0D0/dsc(i)
896 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
897 censc(1,1,-i)=censc(1,1,i)
898 censc(2,1,-i)=censc(2,1,i)
899 censc(3,1,-i)=-censc(3,1,i)
901 read (irotam,*) bsc(j,i)
902 read (irotam,*) (censc(k,j,i),k=1,3),&
903 ((blower(k,l,j),l=1,k),k=1,3)
904 censc(1,j,-i)=censc(1,j,i)
905 censc(2,j,-i)=censc(2,j,i)
906 censc(3,j,-i)=-censc(3,j,i)
907 ! BSC is amplitude of Gaussian
914 akl=akl+blower(k,m,j)*blower(l,m,j)
918 if (((k.eq.3).and.(l.ne.3)) &
919 .or.((l.eq.3).and.(k.ne.3))) then
920 gaussc(k,l,j,-i)=-akl
921 gaussc(l,k,j,-i)=-akl
933 write (iout,'(/a)') 'Parameters of side-chain local geometry'
937 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
938 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
939 ! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
940 ! write (iout,'(a,f10.4,4(16x,f10.4))')
941 ! & 'Center ',(bsc(j,i),j=1,nlobi)
942 ! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
943 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
944 'log h',(bsc(j,i),j=1,nlobi)
945 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
946 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
953 ! blower(k,l,j)=gaussc(ind,j,i)
958 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
959 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
966 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
967 ! added by Urszula Kozlowska 07/11/2007
969 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
977 read(irotam,*) sc_parmin(j,i)
985 ! Read torsional parameters in old format
987 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
989 read (itorp,*) ntortyp,nterm_old
990 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
991 read (itorp,*) (itortyp(i),i=1,ntyp)
993 !el from energy module--------
994 allocate(v1(nterm_old,ntortyp,ntortyp))
995 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
996 !el---------------------------
1002 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1008 write (iout,'(/a/)') 'Torsional constants:'
1011 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1012 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1020 ! Read torsional parameters
1022 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1024 read (itorp,*) ntortyp
1025 read (itorp,*) (itortyp(i),i=1,ntyp)
1026 write (iout,*) 'ntortyp',ntortyp
1028 !el from energy module---------
1029 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1030 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1032 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1033 allocate(vlor2(maxlor,ntortyp,ntortyp))
1034 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1035 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1037 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1038 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1039 !el---------------------------
1041 do i=-ntortyp,ntortyp
1042 do j=-ntortyp,ntortyp
1048 !el---------------------------
1052 itortyp(i)=-itortyp(-i)
1054 ! write (iout,*) 'ntortyp',ntortyp
1056 do j=-ntortyp+1,ntortyp-1
1057 read (itorp,*) nterm(i,j,iblock),&
1059 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1060 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1063 do k=1,nterm(i,j,iblock)
1064 read (itorp,*) kk,v1(k,i,j,iblock),&
1066 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1067 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1068 v0ij=v0ij+si*v1(k,i,j,iblock)
1071 do k=1,nlor(i,j,iblock)
1072 read (itorp,*) kk,vlor1(k,i,j),&
1073 vlor2(k,i,j),vlor3(k,i,j)
1074 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1077 v0(-i,-j,iblock)=v0ij
1084 write (iout,'(/a/)') 'Torsional constants:'
1087 write (iout,*) 'ityp',i,' jtyp',j
1088 write (iout,*) 'Fourier constants'
1089 do k=1,nterm(i,j,iblock)
1090 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1093 write (iout,*) 'Lorenz constants'
1094 do k=1,nlor(i,j,iblock)
1095 write (iout,'(3(1pe15.5))') &
1096 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1103 ! 6/23/01 Read parameters for double torsionals
1105 !el from energy module------------
1106 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1107 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1108 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1109 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1110 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1111 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1112 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1113 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1114 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1115 !---------------------------------
1119 do j=-ntortyp+1,ntortyp-1
1120 do k=-ntortyp+1,ntortyp-1
1121 read (itordp,'(3a1)') t1,t2,t3
1122 ! write (iout,*) "OK onelett",
1125 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1126 .or. t3.ne.toronelet(k)) then
1127 write (iout,*) "Error in double torsional parameter file",&
1130 call MPI_Finalize(Ierror)
1132 stop "Error in double torsional parameter file"
1134 read (itordp,*) ntermd_1(i,j,k,iblock),&
1135 ntermd_2(i,j,k,iblock)
1136 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1137 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1138 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1139 ntermd_1(i,j,k,iblock))
1140 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1141 ntermd_1(i,j,k,iblock))
1142 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1143 ntermd_1(i,j,k,iblock))
1144 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1145 ntermd_1(i,j,k,iblock))
1146 ! Martix of D parameters for one dimesional foureir series
1147 do l=1,ntermd_1(i,j,k,iblock)
1148 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1149 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1150 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1151 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1152 ! write(iout,*) "whcodze" ,
1153 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1155 read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1156 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1157 v2s(m,l,i,j,k,iblock),&
1158 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1159 ! Martix of D parameters for two dimesional fourier series
1160 do l=1,ntermd_2(i,j,k,iblock)
1162 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1163 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1164 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1165 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1174 write (iout,*) 'Constants for double torsionals'
1177 do j=-ntortyp+1,ntortyp-1
1178 do k=-ntortyp+1,ntortyp-1
1179 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1180 ' nsingle',ntermd_1(i,j,k,iblock),&
1181 ' ndouble',ntermd_2(i,j,k,iblock)
1183 write (iout,*) 'Single angles:'
1184 do l=1,ntermd_1(i,j,k,iblock)
1185 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1186 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1187 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1188 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1191 write (iout,*) 'Pairs of angles:'
1192 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1193 do l=1,ntermd_2(i,j,k,iblock)
1194 write (iout,'(i5,20f10.5)') &
1195 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1198 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1199 do l=1,ntermd_2(i,j,k,iblock)
1200 write (iout,'(i5,20f10.5)') &
1201 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1202 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1211 !elwrite(iout,*) "parmread kontrol sc-bb"
1212 ! Read of Side-chain backbone correlation parameters
1213 ! Modified 11 May 2012 by Adasko
1216 read (isccor,*) nsccortyp
1219 !c maxinter is maximum interaction sites
1220 !write(iout,*)"maxterm_sccor",maxterm_sccor
1221 !el from module energy-------------
1222 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1223 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1224 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1225 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1226 !-----------------------------------
1227 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1228 !-----------------------------------
1229 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1230 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1231 -nsccortyp:nsccortyp))
1232 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1233 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1234 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1235 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1236 !-----------------------------------
1237 do i=-nsccortyp,nsccortyp
1238 do j=-nsccortyp,nsccortyp
1242 !-----------------------------------
1244 read (isccor,*) (isccortyp(i),i=1,ntyp)
1246 isccortyp(i)=-isccortyp(-i)
1248 iscprol=isccortyp(20)
1249 ! write (iout,*) 'ntortyp',ntortyp
1251 !c maxinter is maximum interaction sites
1256 nterm_sccor(i,j),nlor_sccor(i,j)
1262 nterm_sccor(-i,j)=nterm_sccor(i,j)
1263 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1264 nterm_sccor(i,-j)=nterm_sccor(i,j)
1265 do k=1,nterm_sccor(i,j)
1266 read (isccor,*) kk,v1sccor(k,l,i,j),&
1268 if (j.eq.iscprol) then
1269 if (i.eq.isccortyp(10)) then
1270 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1271 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1273 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1274 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1275 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1276 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1277 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1278 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1279 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1280 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1283 if (i.eq.isccortyp(10)) then
1284 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1285 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1287 if (j.eq.isccortyp(10)) then
1288 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1289 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1291 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1292 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1293 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1294 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1295 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1296 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1300 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1301 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1302 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1303 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1306 do k=1,nlor_sccor(i,j)
1307 read (isccor,*) kk,vlor1sccor(k,i,j),&
1308 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1309 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1310 (1+vlor3sccor(k,i,j)**2)
1312 v0sccor(l,i,j)=v0ijsccor
1313 v0sccor(l,-i,j)=v0ijsccor1
1314 v0sccor(l,i,-j)=v0ijsccor2
1315 v0sccor(l,-i,-j)=v0ijsccor3
1321 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1324 write (iout,*) 'ityp',i,' jtyp',j
1325 write (iout,*) 'Fourier constants'
1326 do k=1,nterm_sccor(i,j)
1327 write (iout,'(2(1pe15.5))') &
1328 (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
1330 write (iout,*) 'Lorenz constants'
1331 do k=1,nlor_sccor(i,j)
1332 write (iout,'(3(1pe15.5))') &
1333 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1339 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1340 ! interaction energy of the Gly, Ala, and Pro prototypes.
1342 read (ifourier,*) nloctyp
1343 !el write(iout,*)"nloctyp",nloctyp
1344 !el from module energy-------
1345 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1346 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1347 allocate(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1348 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1349 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1350 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1351 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1352 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1354 do ii=-nloctyp-1,nloctyp+1
1362 ctilde(j,i,ii)=0.0d0
1363 dtilde(j,i,ii)=0.0d0
1367 !--------------------------------
1368 allocate(b(13,0:nloctyp))
1372 read (ifourier,*) (b(ii,i),ii=1,13)
1374 write (iout,*) 'Type',i
1375 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1383 B1tilde(1,i) = b(3,i)
1384 B1tilde(2,i) =-b(5,i)
1385 B1tilde(1,-i) =-b(3,i)
1386 B1tilde(2,-i) =b(5,i)
1387 ! b1tilde(1,i)=0.0d0
1388 ! b1tilde(2,i)=0.0d0
1408 Ctilde(1,1,i)=b(7,i)
1409 Ctilde(1,2,i)=b(9,i)
1410 Ctilde(2,1,i)=-b(9,i)
1411 Ctilde(2,2,i)=b(7,i)
1412 Ctilde(1,1,-i)=b(7,i)
1413 Ctilde(1,2,-i)=-b(9,i)
1414 Ctilde(2,1,-i)=b(9,i)
1415 Ctilde(2,2,-i)=b(7,i)
1417 ! Ctilde(1,1,i)=0.0d0
1418 ! Ctilde(1,2,i)=0.0d0
1419 ! Ctilde(2,1,i)=0.0d0
1420 ! Ctilde(2,2,i)=0.0d0
1433 Dtilde(1,1,i)=b(6,i)
1434 Dtilde(1,2,i)=b(8,i)
1435 Dtilde(2,1,i)=-b(8,i)
1436 Dtilde(2,2,i)=b(6,i)
1437 Dtilde(1,1,-i)=b(6,i)
1438 Dtilde(1,2,-i)=-b(8,i)
1439 Dtilde(2,1,-i)=b(8,i)
1440 Dtilde(2,2,-i)=b(6,i)
1442 ! Dtilde(1,1,i)=0.0d0
1443 ! Dtilde(1,2,i)=0.0d0
1444 ! Dtilde(2,1,i)=0.0d0
1445 ! Dtilde(2,2,i)=0.0d0
1446 EE(1,1,i)= b(10,i)+b(11,i)
1447 EE(2,2,i)=-b(10,i)+b(11,i)
1448 EE(2,1,i)= b(12,i)-b(13,i)
1449 EE(1,2,i)= b(12,i)+b(13,i)
1450 EE(1,1,-i)= b(10,i)+b(11,i)
1451 EE(2,2,-i)=-b(10,i)+b(11,i)
1452 EE(2,1,-i)=-b(12,i)+b(13,i)
1453 EE(1,2,-i)=-b(12,i)-b(13,i)
1459 ! ee(2,1,i)=ee(1,2,i)
1464 write (iout,*) 'Type',i
1466 ! write (iout,'(f10.5)') B1(:,i)
1467 write(iout,*) B1(1,i),B1(2,i)
1469 ! write (iout,'(f10.5)') B2(:,i)
1470 write(iout,*) B2(1,i),B2(2,i)
1473 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1477 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1481 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1486 ! Read electrostatic-interaction parameters
1489 write (iout,'(/a)') 'Electrostatic interaction constants:'
1490 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1491 'IT','JT','APP','BPP','AEL6','AEL3'
1493 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1494 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1495 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1496 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1501 app (i,j)=epp(i,j)*rri*rri
1502 bpp (i,j)=-2.0D0*epp(i,j)*rri
1503 ael6(i,j)=elpp6(i,j)*4.2D0**6
1504 ael3(i,j)=elpp3(i,j)*4.2D0**3
1505 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1510 ! Read side-chain interaction parameters.
1512 !el from module energy - COMMON.INTERACT-------
1513 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1514 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1515 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1516 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1517 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1528 !--------------------------------
1530 read (isidep,*) ipot,expon
1531 !el if (ipot.lt.1 .or. ipot.gt.5) then
1532 ! write (iout,'(2a)') 'Error while reading SC interaction',&
1533 ! 'potential file - unknown potential type.'
1537 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1538 ', exponents are ',expon,2*expon
1539 ! goto (10,20,30,30,40) ipot
1541 !----------------------- LJ potential ---------------------------------
1543 ! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1544 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1546 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1547 write (iout,'(a/)') 'The epsilon array:'
1548 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1549 write (iout,'(/a)') 'One-body parameters:'
1550 write (iout,'(a,4x,a)') 'residue','sigma'
1551 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1554 !----------------------- LJK potential --------------------------------
1556 ! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1557 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1558 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1560 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1561 write (iout,'(a/)') 'The epsilon array:'
1562 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1563 write (iout,'(/a)') 'One-body parameters:'
1564 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1565 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1569 !---------------------- GB or BP potential -----------------------------
1573 read (isidep,*)(eps(i,j),j=i,ntyp)
1575 read (isidep,*)(sigma0(i),i=1,ntyp)
1576 read (isidep,*)(sigii(i),i=1,ntyp)
1577 read (isidep,*)(chip(i),i=1,ntyp)
1578 read (isidep,*)(alp(i),i=1,ntyp)
1579 ! For the GB potential convert sigma'**2 into chi'
1582 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1586 write (iout,'(/a/)') 'Parameters of the BP potential:'
1587 write (iout,'(a/)') 'The epsilon array:'
1588 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1589 write (iout,'(/a)') 'One-body parameters:'
1590 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
1592 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
1593 chip(i),alp(i),i=1,ntyp)
1596 !--------------------- GBV potential -----------------------------------
1598 ! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1599 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1600 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
1601 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1603 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1604 write (iout,'(a/)') 'The epsilon array:'
1605 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1606 write (iout,'(/a)') 'One-body parameters:'
1607 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
1608 's||/s_|_^2',' chip ',' alph '
1609 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
1610 sigii(i),chip(i),alp(i),i=1,ntyp)
1613 write (iout,'(2a)') 'Error while reading SC interaction',&
1614 'potential file - unknown potential type.'
1620 !-----------------------------------------------------------------------
1621 ! Calculate the "working" parameters of SC interactions.
1623 !el from module energy - COMMON.INTERACT-------
1624 allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
1625 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
1635 !--------------------------------
1644 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1645 sigma(j,i)=sigma(i,j)
1646 rs0(i,j)=dwa16*sigma(i,j)
1650 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
1651 'Working parameters of the SC interactions:',&
1652 ' a ',' b ',' augm ',' sigma ',' r0 ',&
1657 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1666 sigeps=dsign(1.0D0,epsij)
1668 aa(i,j)=epsij*rrij*rrij
1669 bb(i,j)=-sigeps*epsij*rrij
1673 sigt1sq=sigma0(i)**2
1674 sigt2sq=sigma0(j)**2
1677 ratsig1=sigt2sq/sigt1sq
1678 ratsig2=1.0D0/ratsig1
1679 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1680 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1681 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1685 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1686 sigmaii(i,j)=rsum_max
1687 sigmaii(j,i)=rsum_max
1689 ! sigmaii(i,j)=r0(i,j)
1690 ! sigmaii(j,i)=r0(i,j)
1692 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1693 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1694 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1695 augm(i,j)=epsij*r_augm**(2*expon)
1696 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1703 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
1704 restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
1705 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1710 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
1718 ! Define the SC-p interaction constants
1728 !elwrite(iout,*) "parmread kontrol before oldscp"
1730 ! Define the SC-p interaction constants
1734 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1736 ! aad(i,1)=0.3D0*4.0D0**12
1737 ! Following line for constants currently implemented
1738 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
1739 aad(i,1)=1.5D0*4.0D0**12
1740 ! aad(i,1)=0.17D0*5.6D0**12
1742 ! "Soft" SC-p repulsion
1744 ! Following line for constants currently implemented
1745 ! aad(i,1)=0.3D0*4.0D0**6
1746 ! "Hard" SC-p repulsion
1747 bad(i,1)=3.0D0*4.0D0**6
1748 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1757 ! 8/9/01 Read the SC-p interaction constants from file
1760 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1763 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1764 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1765 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1766 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1770 write (iout,*) "Parameters of SC-p interactions:"
1772 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
1773 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1778 ! Define the constants of the disulfide bridge
1782 ! Old arbitrary potential - commented out.
1787 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1788 ! energy surface of diethyl disulfide.
1789 ! A. Liwo and U. Kozlowska, 11/24/03
1800 write (iout,'(/a)') "Disulfide bridge parameters:"
1801 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1802 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1803 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1804 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
1808 end subroutine parmread
1810 !-----------------------------------------------------------------------------
1812 !-----------------------------------------------------------------------------
1813 subroutine mygetenv(string,var)
1817 ! This subroutine passes the environmental variables to FORTRAN program.
1818 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
1819 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
1820 ! reads the environmental variables from $HOME/.env
1822 ! Usage: As for the standard FORTRAN GETENV subroutine.
1824 ! Purpose: some versions/installations of MPI do not transfer the environmental
1825 ! variables to slave processors, if these variables are set in the shell script
1826 ! from which mpirun is called.
1835 character*(*) :: string,var
1836 #if defined(MYGETENV) && defined(MPI)
1837 ! include "DIMENSIONS.ZSCOPT"
1839 ! include "COMMON.MPI"
1840 !el character*360 ucase
1842 character(len=360) :: string1(360),karta
1843 character(len=240) :: home
1846 call getenv("HOME",home)
1847 open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
1849 read (99,end=111,err=111,'(a)') karta
1853 call split_string(karta,string1,80,n)
1854 if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
1855 string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
1857 print *,"Processor",me,": ",var(:ilen(var)),&
1858 " assigned to ",string(:ilen(string))
1863 111 print *,"Environment variable ",string(:ilen(string))," not set."
1866 112 print *,"Error opening environment file!"
1868 call getenv(string,var)
1871 end subroutine mygetenv
1872 !-----------------------------------------------------------------------------
1874 !-----------------------------------------------------------------------------
1875 subroutine read_general_data(*)
1877 use control_data, only:indpdb,symetr
1878 use energy_data, only:distchainmax
1880 ! include "DIMENSIONS"
1881 ! include "DIMENSIONS.ZSCOPT"
1882 ! include "DIMENSIONS.FREE"
1883 ! include "COMMON.TORSION"
1884 ! include "COMMON.INTERACT"
1885 ! include "COMMON.IOUNITS"
1886 ! include "COMMON.TIME1"
1887 ! include "COMMON.PROT"
1888 ! include "COMMON.PROTFILES"
1889 ! include "COMMON.CHAIN"
1890 ! include "COMMON.NAMES"
1891 ! include "COMMON.FFIELD"
1892 ! include "COMMON.ENEPS"
1893 ! include "COMMON.WEIGHTS"
1894 ! include "COMMON.FREE"
1895 ! include "COMMON.CONTROL"
1896 ! include "COMMON.ENERGIES"
1897 character(len=800) :: controlcard
1898 integer :: i,j,k,ii,n_ene_found
1899 integer :: ind,itype1,itype2,itypf,itypsc,itypp
1902 !el character*16 ucase
1903 character(len=16) :: key
1905 call card_concat(controlcard,.true.)
1906 call readi(controlcard,"N_ENE",n_eneW,max_eneW)
1907 if (n_eneW.gt.max_eneW) then
1908 write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
1912 call readi(controlcard,"NPARMSET",nparmset,1)
1913 !elwrite(iout,*)"in read_gen data"
1914 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
1915 call readi(controlcard,"IPARMPRINT",iparmprint,1)
1916 write (iout,*) "PARMPRINT",iparmprint
1917 if (nparmset.gt.max_parm) then
1918 write (iout,*) "Error: parameter out of range: NPARMSET",&
1922 !elwrite(iout,*)"in read_gen data"
1923 call readi(controlcard,"MAXIT",maxit,5000)
1924 call reada(controlcard,"FIMIN",fimin,1.0d-3)
1925 call readi(controlcard,"ENSEMBLES",ensembles,0)
1926 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
1927 write (iout,*) "Number of energy parameter sets",nparmset
1928 allocate(isampl(nparmset))
1929 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
1930 write (iout,*) "MaxSlice",MaxSlice
1931 call readi(controlcard,"NSLICE",nslice,1)
1932 !elwrite(iout,*)"in read_gen data"
1934 if (nslice.gt.MaxSlice) then
1935 write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
1939 write (iout,*) "Frequency of storing conformations",&
1940 (isampl(i),i=1,nparmset)
1941 write (iout,*) "Maxit",maxit," Fimin",fimin
1942 call readi(controlcard,"NQ",nQ,1)
1943 if (nQ.gt.MaxQ) then
1944 write (iout,*) "Error: parameter out of range: NQ",nq,&
1949 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
1950 call reada(controlcard,"DELTA",delta,1.0d-2)
1951 call readi(controlcard,"EINICHECK",einicheck,2)
1952 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
1953 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
1954 call readi(controlcard,"RESCALE",rescale_modeW,1)
1955 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
1956 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1957 call readi(controlcard,'SYM',symetr,1)
1958 write (iout,*) "DISTCHAINMAX",distchainmax
1959 write (iout,*) "delta",delta
1960 write (iout,*) "einicheck",einicheck
1961 write (iout,*) "rescale_mode",rescale_modeW
1963 bxfile=index(controlcard,"BXFILE").gt.0
1964 cxfile=index(controlcard,"CXFILE").gt.0
1965 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
1967 histfile=index(controlcard,"HISTFILE").gt.0
1968 histout=index(controlcard,"HISTOUT").gt.0
1969 entfile=index(controlcard,"ENTFILE").gt.0
1970 zscfile=index(controlcard,"ZSCFILE").gt.0
1971 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1972 write (iout,*) "with_dihed_constr ",with_dihed_constr
1973 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
1975 end subroutine read_general_data
1976 !------------------------------------------------------------------------------
1977 subroutine read_efree(*)
1979 ! Read molecular data
1982 ! include 'DIMENSIONS'
1983 ! include 'DIMENSIONS.ZSCOPT'
1984 ! include 'DIMENSIONS.COMPAR'
1985 ! include 'DIMENSIONS.FREE'
1986 ! include 'COMMON.IOUNITS'
1987 ! include 'COMMON.TIME1'
1988 ! include 'COMMON.SBRIDGE'
1989 ! include 'COMMON.CONTROL'
1990 ! include 'COMMON.CHAIN'
1991 ! include 'COMMON.HEADER'
1992 ! include 'COMMON.GEO'
1993 ! include 'COMMON.FREE'
1994 character(len=320) :: controlcard !,ucase
1995 integer :: iparm,ib,i,j,npars
2005 ! call alloc_wham_arrays
2006 ! allocate(nT_h(nParmSet))
2007 ! allocate(replica(nParmSet))
2008 ! allocate(umbrella(nParmSet))
2009 ! allocate(read_iset(nParmSet))
2010 ! allocate(nT_h(nParmSet))
2014 call card_concat(controlcard,.true.)
2015 call readi(controlcard,'NT',nT_h(iparm),1)
2016 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
2018 if (nT_h(iparm).gt.MaxT_h) then
2019 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
2023 replica(iparm)=index(controlcard,"REPLICA").gt.0
2024 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
2025 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
2026 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
2027 replica(iparm)," umbrella ",umbrella(iparm),&
2028 " read_iset",read_iset(iparm)
2031 call card_concat(controlcard,.true.)
2032 call readi(controlcard,'NR',nR(ib,iparm),1)
2033 if (umbrella(iparm)) then
2036 nRR(ib,iparm)=nR(ib,iparm)
2038 if (nR(ib,iparm).gt.MaxR) then
2039 write (iout,*) "Error: parameter out of range: NR",&
2043 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
2044 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
2045 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
2048 call card_concat(controlcard,.true.)
2049 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
2051 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
2056 write (iout,*) "ib",ib," beta_h",&
2057 1.0d0/(0.001987*beta_h(ib,iparm))
2058 write (iout,*) "nR",nR(ib,iparm)
2059 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
2061 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
2062 "q0",(q0(j,i,ib,iparm),j=1,nQ)
2074 nR(ib,iparm)=nR(ib,1)
2075 if (umbrella(iparm)) then
2078 nRR(ib,iparm)=nR(ib,1)
2080 beta_h(ib,iparm)=beta_h(ib,1)
2082 f(i,ib,iparm)=f(i,ib,1)
2084 KH(j,i,ib,iparm)=KH(j,i,ib,1)
2085 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
2088 replica(iparm)=replica(1)
2089 umbrella(iparm)=umbrella(1)
2090 read_iset(iparm)=read_iset(1)
2097 end subroutine read_efree
2098 !-----------------------------------------------------------------------------
2099 subroutine read_protein_data(*)
2101 ! include "DIMENSIONS"
2102 ! include "DIMENSIONS.ZSCOPT"
2103 ! include "DIMENSIONS.FREE"
2107 integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
2108 ! include "COMMON.MPI"
2110 ! include "COMMON.CHAIN"
2111 ! include "COMMON.IOUNITS"
2112 ! include "COMMON.PROT"
2113 ! include "COMMON.PROTFILES"
2114 ! include "COMMON.NAMES"
2115 ! include "COMMON.FREE"
2116 ! include "COMMON.OBCINKA"
2117 character(len=64) :: nazwa
2118 character(len=16000) :: controlcard
2119 integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
2120 !el external ilen,iroof
2129 ! Read names of files with conformation data.
2130 if (replica(iparm)) then
2136 do ii=1,nRR(ib,iparm)
2137 write (iout,*) "Parameter set",iparm," temperature",ib,&
2140 call card_concat(controlcard,.true.)
2141 write (iout,*) controlcard(:ilen(controlcard))
2142 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
2143 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
2144 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
2145 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
2146 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
2147 maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
2148 call reada(controlcard,"TIME_START",&
2149 time_start_collect(ii,ib,iparm),0.0d0)
2150 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
2152 write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
2153 " rec_end",rec_end(ii,ib,iparm)
2154 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
2155 " time_end",time_end_collect(ii,ib,iparm)
2157 if (replica(iparm)) then
2158 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
2159 write (iout,*) "Number of trajectories",totraj(ii,iparm)
2162 if (nfile_bin(ii,ib,iparm).lt.2 &
2163 .and. nfile_asc(ii,ib,iparm).eq.0 &
2164 .and. nfile_cx(ii,ib,iparm).eq.0) then
2165 write (iout,*) "Error - no action specified!"
2168 if (nfile_bin(ii,ib,iparm).gt.0) then
2169 call card_concat(controlcard,.false.)
2170 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
2171 maxfile_prot,nfile_bin(ii,ib,iparm))
2173 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
2174 write(iout,*) (protfiles(i,1,ii,ib,iparm),&
2175 i=1,nfile_bin(ii,ib,iparm))
2178 if (nfile_asc(ii,ib,iparm).gt.0) then
2179 call card_concat(controlcard,.false.)
2180 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2181 maxfile_prot,nfile_asc(ii,ib,iparm))
2183 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
2184 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2185 i=1,nfile_asc(ii,ib,iparm))
2187 else if (nfile_cx(ii,ib,iparm).gt.0) then
2188 call card_concat(controlcard,.false.)
2189 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2190 maxfile_prot,nfile_cx(ii,ib,iparm))
2192 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
2193 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2194 i=1,nfile_cx(ii,ib,iparm))
2203 end subroutine read_protein_data
2204 !-------------------------------------------------------------------------------
2205 subroutine readsss(rekord,lancuch,wartosc,default)
2207 character*(*) :: rekord,lancuch,wartosc,default
2208 character(len=80) :: aux
2209 integer :: lenlan,lenrec,iread,ireade
2213 lenlan=ilen(lancuch)
2215 iread=index(rekord,lancuch(:lenlan)//"=")
2216 ! print *,"rekord",rekord," lancuch",lancuch
2217 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2218 if (iread.eq.0) then
2222 iread=iread+lenlan+1
2223 ! print *,"iread",iread
2224 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2225 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2227 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2229 ! print *,"iread",iread
2230 if (iread.gt.lenrec) then
2235 ! print *,"ireade",ireade
2236 do while (ireade.lt.lenrec .and. &
2237 .not.iblnk(rekord(ireade:ireade)))
2240 wartosc=rekord(iread:ireade)
2242 end subroutine readsss
2243 !----------------------------------------------------------------------------
2244 subroutine multreads(rekord,lancuch,tablica,dim,default)
2247 character*(*) rekord,lancuch,tablica(dim),default
2248 character(len=80) :: aux
2249 integer :: lenlan,lenrec,iread,ireade
2256 lenlan=ilen(lancuch)
2258 iread=index(rekord,lancuch(:lenlan)//"=")
2259 ! print *,"rekord",rekord," lancuch",lancuch
2260 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2261 if (iread.eq.0) return
2262 iread=iread+lenlan+1
2264 ! print *,"iread",iread
2265 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2266 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2268 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2270 ! print *,"iread",iread
2271 if (iread.gt.lenrec) return
2273 ! print *,"ireade",ireade
2274 do while (ireade.lt.lenrec .and. &
2275 .not.iblnk(rekord(ireade:ireade)))
2278 tablica(i)=rekord(iread:ireade)
2281 end subroutine multreads
2282 !----------------------------------------------------------------------------
2283 subroutine split_string(rekord,tablica,dim,nsub)
2285 integer :: dim,nsub,i,ii,ll,kk
2286 character*(*) tablica(dim)
2287 character*(*) rekord
2297 ! Find the start of term name
2299 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2302 ! Parse the name into TABLICA(i) until blank found
2303 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2305 tablica(i)(kk:kk)=rekord(ii:ii)
2308 if (kk.gt.0) nsub=nsub+1
2309 if (ii.gt.ll) return
2312 end subroutine split_string
2313 !--------------------------------------------------------------------------------
2315 !--------------------------------------------------------------------------------
2316 subroutine read_compar
2318 ! Read molecular data
2320 use conform_compar, only:alloc_compar_arrays
2321 use control_data, only:pdbref
2322 use geometry_data, only:deg2rad,rad2deg
2324 ! include 'DIMENSIONS'
2325 ! include 'DIMENSIONS.ZSCOPT'
2326 ! include 'DIMENSIONS.COMPAR'
2327 ! include 'DIMENSIONS.FREE'
2328 ! include 'COMMON.IOUNITS'
2329 ! include 'COMMON.TIME1'
2330 ! include 'COMMON.SBRIDGE'
2331 ! include 'COMMON.CONTROL'
2332 ! include 'COMMON.COMPAR'
2333 ! include 'COMMON.CHAIN'
2334 ! include 'COMMON.HEADER'
2335 ! include 'COMMON.GEO'
2336 ! include 'COMMON.FREE'
2337 character(len=320) :: controlcard !,ucase
2338 character(len=64) :: wfile
2342 !elwrite(iout,*)"jestesmy w read_compar"
2343 call card_concat(controlcard,.true.)
2344 pdbref=(index(controlcard,'PDBREF').gt.0)
2345 call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
2346 call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
2347 call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
2348 call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
2349 verbose = index(controlcard,"VERBOSE").gt.0
2350 lgrp=index(controlcard,"STATIN").gt.0
2351 lgrp_out=index(controlcard,"STATOUT").gt.0
2352 merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
2353 binary = index(controlcard,"BINARY").gt.0
2354 rmscut_base_up=rmscut_base_up/50
2355 rmscut_base_low=rmscut_base_low/50
2356 call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
2357 call readi(controlcard,'NLEVEL',nlevel,1)
2358 if (nlevel.lt.0) then
2360 call alloc_compar_arrays(maxfrag,1)
2363 allocate(nfrag(nlevel))
2365 ! Read the data pertaining to elementary fragments (level 1)
2366 call readi(controlcard,'NFRAG',nfrag(1),0)
2367 write(iout,*)"nfrag(1)",nfrag(1)
2368 call alloc_compar_arrays(nfrag(1),nlevel)
2370 call card_concat(controlcard,.true.)
2371 write (iout,*) controlcard(:ilen(controlcard))
2372 call readi(controlcard,'NPIECE',npiece(j,1),0)
2373 call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
2374 call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
2375 call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
2376 call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
2377 call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
2378 call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
2379 call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
2380 call readi(controlcard,'RMS',irms(j,1),0)
2381 call readi(controlcard,'LOCAL',iloc(j),1)
2382 call readi(controlcard,'ELCONT',ielecont(j,1),1)
2383 if (ielecont(j,1).eq.0) then
2384 call readi(controlcard,'SCCONT',isccont(j,1),1)
2386 ang_cut(j)=ang_cut(j)*deg2rad
2387 ang_cut1(j)=ang_cut1(j)*deg2rad
2389 call card_concat(controlcard,.true.)
2390 call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
2391 call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
2393 write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
2394 (ifrag(1,k,j),ifrag(2,k,j),&
2395 k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
2396 " ang_cut1",ang_cut1(j)*rad2deg
2397 write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
2398 write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
2399 write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
2400 " ilocal",iloc(j)," isccont",isccont(j,1)
2402 ! Read data pertaning to higher levels
2404 call card_concat(controlcard,.true.)
2405 call readi(controlcard,'NFRAG',NFRAG(i),0)
2406 write (iout,*) "i",i," nfrag",nfrag(i)
2408 call card_concat(controlcard,.true.)
2410 call readi(controlcard,'ELCONT',ielecont(j,i),0)
2411 if (ielecont(j,i).eq.0) then
2412 call readi(controlcard,'SCCONT',isccont(j,i),1)
2414 call readi(controlcard,'RMS',irms(j,i),0)
2420 call readi(controlcard,'NPIECE',npiece(j,i),0)
2421 call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
2422 call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
2423 call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
2425 call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
2426 call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
2427 write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
2428 n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
2429 " isccont",isccont(j,i)," irms",irms(j,i)
2430 write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
2431 write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
2432 write(iout,*)"nc_frac",nc_fragm(j,i),&
2433 " nc_req",nc_req_setf(j,i)
2436 if (binary) write (iout,*) "Classes written in binary format."
2439 call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
2440 call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
2441 call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
2442 call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
2443 call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
2444 call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
2445 call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
2446 call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
2447 call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
2448 call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
2449 call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
2450 call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
2451 call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
2452 call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
2453 call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
2454 call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
2455 call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
2456 call readi(controlcard,'RMS_SINGLE',irms_single,0)
2457 call readi(controlcard,'CONT_SINGLE',icont_single,1)
2458 call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
2459 call readi(controlcard,'RMS_PAIR',irms_pair,0)
2460 call readi(controlcard,'CONT_PAIR',icont_pair,1)
2461 call readi(controlcard,'SPLIT_BET',isplit_bet,0)
2462 angcut_hel=angcut_hel*deg2rad
2463 angcut1_hel=angcut1_hel*deg2rad
2464 angcut_bet=angcut_bet*deg2rad
2465 angcut1_bet=angcut1_bet*deg2rad
2466 angcut_strand=angcut_strand*deg2rad
2467 angcut1_strand=angcut1_strand*deg2rad
2468 write (iout,*) "Automatic detection of structural elements"
2469 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2470 ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2471 ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
2472 ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2473 ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
2474 ' SPLIT_BET',isplit_bet
2475 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2476 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2477 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2478 ' MAXANG_HEL',angcut1_hel*rad2deg
2479 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2480 ' MAXANG_BET',angcut1_bet*rad2deg
2481 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2482 ' MAXANG_STRAND',angcut1_strand*rad2deg
2483 write (iout,*) 'FRAC_MIN',frac_min_set
2485 end subroutine read_compar
2486 !--------------------------------------------------------------------------------
2488 !--------------------------------------------------------------------------------
2489 subroutine read_ref_structure(*)
2491 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
2494 use control_data, only:pdbref
2495 use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
2496 nstart_sup,nstart_seq,nperm,nres0
2497 use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype
2498 use compare, only:seq_comp !,contact,elecont
2499 use geometry, only:chainbuild,dist
2500 use io_config, only:readpdb
2502 use conform_compar, only:contact,elecont
2504 ! include 'DIMENSIONS'
2505 ! include 'DIMENSIONS.ZSCOPT'
2506 ! include 'DIMENSIONS.COMPAR'
2507 ! include 'COMMON.IOUNITS'
2508 ! include 'COMMON.GEO'
2509 ! include 'COMMON.VAR'
2510 ! include 'COMMON.INTERACT'
2511 ! include 'COMMON.LOCAL'
2512 ! include 'COMMON.NAMES'
2513 ! include 'COMMON.CHAIN'
2514 ! include 'COMMON.FFIELD'
2515 ! include 'COMMON.SBRIDGE'
2516 ! include 'COMMON.HEADER'
2517 ! include 'COMMON.CONTROL'
2518 ! include 'COMMON.CONTACTS1'
2519 ! include 'COMMON.PEPTCONT'
2520 ! include 'COMMON.TIME1'
2521 ! include 'COMMON.COMPAR'
2522 character(len=4) :: sequence(nres)
2524 !el real(kind=8) :: x(maxvar)
2525 integer :: itype_pdb(nres)
2526 !el logical seq_comp
2527 integer :: i,j,k,nres_pdb,iaux
2528 real(kind=8) :: ddsc !el,dist
2529 integer :: kkk !,ilen
2533 write (iout,*) "pdbref",pdbref
2535 read(inp,'(a)') pdbfile
2536 write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
2537 pdbfile(:ilen(pdbfile))
2538 open(ipdbin,file=pdbfile,status='old',err=33)
2540 33 write (iout,'(a)') 'Error opening PDB file.'
2544 itype_pdb(i)=itype(i)
2551 itype_pdb(i)=itype(i)
2559 if (nsup.le.(nct-nnt+1)) then
2560 do i=0,nct-nnt+1-nsup
2561 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),&
2563 do j=nnt+nsup-1,nnt,-1
2565 cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
2568 do j=nnt+nsup-1,nnt,-1
2570 cref(k,j+i,kkk)=cref(k,j,kkk)
2572 phi_ref(j+i)=phi_ref(j)
2573 theta_ref(j+i)=theta_ref(j)
2574 alph_ref(j+i)=alph_ref(j)
2575 omeg_ref(j+i)=omeg_ref(j)
2579 write (iout,'(i5,3f10.5,5x,3f10.5)') &
2580 j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
2588 write (iout,'(a)') &
2589 'Error - sequences to be superposed do not match.'
2592 do i=0,nsup-(nct-nnt+1)
2593 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),&
2596 nstart_sup=nstart_sup+i
2601 write (iout,'(a)') &
2602 'Error - sequences to be superposed do not match.'
2606 write (iout,'(a,i5)') &
2607 'Experimental structure begins at residue',nstart_seq
2609 call read_angles(inp,*38)
2611 38 write (iout,'(a)') 'Error reading reference structure.'
2620 cref(j,i,kkk)=c(j,i)
2624 nend_sup=nstart_sup+nsup-1
2627 c(j,i)=cref(j,i,kkk)
2632 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
2634 if (itype(i).ne.10) then
2635 ddsc = dist(i,nres+i)
2637 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
2641 dc_norm(j,nres+i)=0.0d0
2644 ! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
2645 ! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
2646 ! dc_norm(3,nres+i)**2
2648 dc(j,i)=c(j,i+1)-c(j,i)
2652 dc_norm(j,i)=dc(j,i)/ddsc
2655 ! print *,"Calling contact"
2656 call contact(.true.,ncont_ref,icont_ref(1,1),&
2657 nstart_sup,nend_sup)
2658 ! print *,"Calling elecont"
2659 call elecont(.true.,ncont_pept_ref,&
2660 icont_pept_ref(1,1),&
2661 nstart_sup,nend_sup)
2662 write (iout,'(a,i3,a,i3,a,i3,a)') &
2663 'Number of residues to be superposed:',nsup,&
2664 ' (from residue',nstart_sup,' to residue',&
2667 end subroutine read_ref_structure
2668 !--------------------------------------------------------------------------------
2670 !--------------------------------------------------------------------------------
2671 subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
2673 use geometry_data, only:nres,c
2674 use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype
2675 ! implicit real*8 (a-h,o-z)
2676 ! include 'DIMENSIONS'
2677 ! include 'DIMENSIONS.ZSCOPT'
2678 ! include 'COMMON.CHAIN'
2679 ! include 'COMMON.INTERACT'
2680 ! include 'COMMON.NAMES'
2681 ! include 'COMMON.IOUNITS'
2682 ! include 'COMMON.HEADER'
2683 ! include 'COMMON.SBRIDGE'
2684 character(len=50) :: tytul
2685 character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
2686 'D','E','F','G','H','I','J'/),shape(chainid))
2687 integer,dimension(nres) :: ica !(maxres)
2688 real(kind=8) :: temp,efree,etot,entropy,rmsdev
2689 integer :: ii,i,j,iti,ires,iatom,ichain
2690 write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
2692 write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
2694 write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
2701 if (iti.eq.ntyp1) then
2704 write (ipdb,'(a)') 'TER'
2709 write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
2713 write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
2714 ires,(c(j,nres+i),j=1,3)
2718 write (ipdb,'(a)') 'TER'
2720 if (itype(i).eq.ntyp1) cycle
2721 if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
2722 write (ipdb,30) ica(i),ica(i+1)
2723 else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
2724 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2725 else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
2726 write (ipdb,30) ica(i),ica(i)+1
2729 if (itype(nct).ne.10) then
2730 write (ipdb,30) ica(nct),ica(nct)+1
2733 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2735 write (ipdb,'(a6)') 'ENDMDL'
2736 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2737 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2738 30 FORMAT ('CONECT',8I5)
2740 end subroutine pdboutW
2742 !------------------------------------------------------------------------------
2744 !-----------------------------------------------------------------------------
2745 !-----------------------------------------------------------------------------