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,nres_molec
135 use control_data, only:iscode
136 use io_base, only:rescode
137 use control, only:setup_var,init_int_table
138 use geometry, only:alloc_geo_arrays
139 use energy, only:alloc_ener_arrays
140 ! implicit real*8 (a-h,o-z)
141 ! include 'DIMENSIONS'
142 ! include 'DIMENSIONS.ZSCOPT'
143 ! include 'COMMON.IOUNITS'
144 ! include 'COMMON.GEO'
145 ! include 'COMMON.VAR'
146 ! include 'COMMON.INTERACT'
147 ! include 'COMMON.LOCAL'
148 ! include 'COMMON.NAMES'
149 ! include 'COMMON.CHAIN'
150 ! include 'COMMON.FFIELD'
151 ! include 'COMMON.SBRIDGE'
152 ! include 'COMMON.TORCNSTR'
153 ! include 'COMMON.CONTROL'
154 character(len=4),dimension(:),allocatable :: sequence !(nres)
155 !el integer :: rescode
156 !el real(kind=8) :: x(maxvar)
157 character(len=320) :: controlcard !,ucase
158 integer,dimension(nres,5) :: itype_pdb !(maxres)
159 integer :: i,j,i1,i2,it1,it2,mnum
160 real(kind=8) :: scalscp
161 !el logical :: seq_comp
162 write(iout,*) "START MOLREAD"
166 print *,"nres_molec, initial",nres_molec(i)
169 call card_concat(controlcard,.true.)
170 call reada(controlcard,'SCAL14',scal14,0.4d0)
171 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
172 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
173 call reada(controlcard,'TEMP0',temp0,300.0d0) !el
174 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
175 r0_corr=cutoff_corr-delt_corr
176 call readi(controlcard,"NRES",nres_molec(1),0)
177 call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
178 call readi(controlcard,"NRES_CAT",nres_molec(5),0)
181 nres=nres_molec(i)+nres
183 ! allocate(sequence(nres+1))
184 !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
185 call alloc_geo_arrays
186 call alloc_ener_arrays
187 ! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
188 !----------------------------
189 allocate(c(3,2*nres+2))
190 allocate(dc(3,0:2*nres+2))
191 allocate(itype(nres+2,5))
192 allocate(itel(nres+2))
193 if (.not. allocated(molnum)) allocate(molnum(nres+2))
209 !--------------------------
211 allocate(sequence(nres+1))
212 iscode=index(controlcard,"ONE_LETTER")
214 write (iout,*) "Error: no residues in molecule"
217 if (nres.gt.maxres) then
218 write (iout,*) "Error: too many residues",nres,maxres
220 write(iout,*) 'nres=',nres
222 ! Read sequence of the protein
223 if (iscode.gt.0) then
224 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
226 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
228 ! Convert sequence to numeric code
232 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
234 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
237 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
239 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
242 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
245 write (iout,*) "Numeric code:"
246 write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres)
250 if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
252 if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
256 else if (iabs(itype(i+1,mnum)).ne.20) then
258 else if (iabs(itype(i,mnum)).ne.20) then
265 write (iout,*) "ITEL"
268 write (iout,*) i,itype(i,mnum),itel(i)
273 if (with_dihed_constr) then
275 read (inp,*) ndih_constr
276 if (ndih_constr.gt.0) then
278 write (iout,*) 'FTORS',ftors
279 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
281 'There are',ndih_constr,' constraints on phi angles.'
283 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
286 phi0(i)=deg2rad*phi0(i)
287 drange(i)=deg2rad*drange(i)
295 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
296 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
297 write(iout,*) 'NNT=',NNT,' NCT=',NCT
301 write (iout,'(/a,i3,a)') 'The chain contains',ns,&
302 ' disulfide-bridging cysteines.'
303 write (iout,'(20i4)') (iss(i),i=1,ns)
304 write (iout,'(/a/)') 'Pre-formed links are:'
311 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
312 restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',&
313 dhpb(i),ebr,forcon(i)
318 end subroutine molread
319 !-----------------------------------------------------------------------------
321 !-----------------------------------------------------------------------------
322 subroutine parmread(iparm,*)
327 ! Read the parameters of the probability distributions of the virtual-bond
328 ! valence angles and the side chains and energy parameters.
334 use control_data, only: maxterm,maxlor,maxterm_sccor,& !maxtor
335 maxtermd_1,maxtermd_2 !,maxthetyp,maxthetyp1
339 use io_config, only: printmat
340 use control, only: getenv_loc
347 ! implicit real*8 (a-h,o-z)
348 ! include 'DIMENSIONS'
349 ! include 'DIMENSIONS.ZSCOPT'
350 ! include 'DIMENSIONS.FREE'
351 ! include 'COMMON.IOUNITS'
352 ! include 'COMMON.CHAIN'
353 ! include 'COMMON.INTERACT'
354 ! include 'COMMON.GEO'
355 ! include 'COMMON.LOCAL'
356 ! include 'COMMON.TORSION'
357 ! include 'COMMON.FFIELD'
358 ! include 'COMMON.NAMES'
359 ! include 'COMMON.SBRIDGE'
360 ! include 'COMMON.WEIGHTS'
361 ! include 'COMMON.ENEPS'
362 ! include 'COMMON.SCCOR'
363 ! include 'COMMON.SCROT'
364 ! include 'COMMON.FREE'
365 character(len=1) :: t1,t2,t3
366 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
367 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
369 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
370 character(len=800) :: controlcard
371 character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
372 tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
376 character(len=16) :: key
378 !el real(kind=8) :: ip,mp
379 real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
380 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
381 real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
383 integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n
384 integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
388 ! Set LPRINT=.TRUE. for debugging
389 dwa16=2.0d0**(1.0d0/6.0d0)
392 ! Assign virtual-bond length
395 vblinv2=vblinv*vblinv
397 call card_concat(controlcard,.true.)
400 allocate(ww(max_eneW))
402 key = wname(i)(:ilen(wname(i)))
403 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
406 write (iout,*) "iparm",iparm," myparm",myparm
407 ! If reading not own parameters, skip assignment
409 if (iparm.eq.myparm .or. .not.separate_parset) then
412 ! Setup weights for UNRES
437 allocate(weights(n_ene))
452 weights(15)=0 !wstrain !
453 weights(16)=0 !wvdwpp !
455 weights(18)=0 !scal14 !
461 call card_concat(controlcard,.false.)
463 ! Return if not own parameters
465 if (iparm.ne.myparm .and. separate_parset) return
467 call reads(controlcard,"BONDPAR",bondname_t,bondname)
468 open (ibond,file=bondname_t,status='old')
470 call reads(controlcard,"THETPAR",thetname_t,thetname)
471 open (ithep,file=thetname_t,status='old')
473 call reads(controlcard,"ROTPAR",rotname_t,rotname)
474 open (irotam,file=rotname_t,status='old')
476 call reads(controlcard,"TORPAR",torname_t,torname)
477 open (itorp,file=torname_t,status='old')
479 call reads(controlcard,"TORDPAR",tordname_t,tordname)
480 open (itordp,file=tordname_t,status='old')
482 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
483 open (isccor,file=sccorname_t,status='old')
485 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
486 open (ifourier,file=fouriername_t,status='old')
488 call reads(controlcard,"ELEPAR",elename_t,elename)
489 open (ielep,file=elename_t,status='old')
491 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
492 open (isidep,file=sidename_t,status='old')
494 call reads(controlcard,"SCPPAR",scpname_t,scpname)
495 open (iscpp,file=scpname_t,status='old')
497 call getenv_loc('IONPAR',ionname)
498 open (iion,file=ionname,status='old')
500 write (iout,*) "Parameter set:",iparm
501 write (iout,*) "Energy-term weights:"
503 write (iout,'(a16,f10.5)') wname(i),ww(i)
505 write (iout,*) "Sidechain potential file : ",&
506 sidename_t(:ilen(sidename_t))
508 write (iout,*) "SCp potential file : ",&
509 scpname_t(:ilen(scpname_t))
511 write (iout,*) "Electrostatic potential file : ",&
512 elename_t(:ilen(elename_t))
513 write (iout,*) "Cumulant coefficient file : ",&
514 fouriername_t(:ilen(fouriername_t))
515 write (iout,*) "Torsional parameter file : ",&
516 torname_t(:ilen(torname_t))
517 write (iout,*) "Double torsional parameter file : ",&
518 tordname_t(:ilen(tordname_t))
519 write (iout,*) "Backbone-rotamer parameter file : ",&
520 sccorname(:ilen(sccorname))
521 write (iout,*) "Bond & inertia constant file : ",&
522 bondname_t(:ilen(bondname_t))
523 write (iout,*) "Bending parameter file : ",&
524 thetname_t(:ilen(thetname_t))
525 write (iout,*) "Rotamer parameter file : ",&
526 rotname_t(:ilen(rotname_t))
529 ! Read the virtual-bond parameters, masses, and moments of inertia
530 ! and Stokes' radii of the peptide group and side chains
532 allocate(dsc(ntyp1)) !(ntyp1)
533 allocate(dsc_inv(ntyp1)) !(ntyp1)
534 allocate(nbondterm(ntyp)) !(ntyp)
535 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
536 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
537 !el allocate(msc(ntyp+1)) !(ntyp+1)
538 !el allocate(isc(ntyp+1)) !(ntyp+1)
539 !el allocate(restok(ntyp+1)) !(ntyp+1)
540 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
543 read (ibond,*) vbldp0,akp
546 read (ibond,*) vbldsc0(1,i),aksc(1,i)
547 dsc(i) = vbldsc0(1,i)
551 dsc_inv(i)=1.0D0/dsc(i)
555 read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
557 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
559 dsc(i) = vbldsc0(1,i)
563 dsc_inv(i)=1.0D0/dsc(i)
568 write(iout,'(/a/)')"Force constants virtual bonds:"
569 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
571 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
573 write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
574 vbldsc0(1,i),aksc(1,i),abond0(1,i)
576 write (iout,'(13x,3f10.5)') &
577 vbldsc0(j,i),aksc(j,i),abond0(j,i)
581 if (.not. allocated(msc)) allocate(msc(ntyp1,5))
582 if (.not. allocated(restok)) allocate(restok(ntyp1,5))
585 read(iion,*) msc(i,5),restok(i,5)
586 print *,msc(i,5),restok(i,5)
590 read (iion,*) ncatprotparm
591 allocate(catprm(ncatprotparm,4))
593 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
598 !----------------------------------------------------
599 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
600 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
601 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
602 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
603 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
604 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
610 athet(j,i,ichir1,ichir2)=0.0D0
611 bthet(j,i,ichir1,ichir2)=0.0D0
625 !elwrite(iout,*) "parmread kontrol"
629 ! Read the parameters of the probability distribution/energy expression
630 ! of the virtual-bond valence angles theta
633 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
634 (bthet(j,i,1,1),j=1,2)
635 read (ithep,*) (polthet(j,i),j=0,3)
636 !elwrite(iout,*) "parmread kontrol in cryst_theta"
637 read (ithep,*) (gthet(j,i),j=1,3)
638 !elwrite(iout,*) "parmread kontrol in cryst_theta"
639 read (ithep,*) theta0(i),sig0(i),sigc0(i)
641 !elwrite(iout,*) "parmread kontrol in cryst_theta"
643 !elwrite(iout,*) "parmread kontrol in cryst_theta"
645 athet(1,i,1,-1)=athet(1,i,1,1)
646 athet(2,i,1,-1)=athet(2,i,1,1)
647 bthet(1,i,1,-1)=-bthet(1,i,1,1)
648 bthet(2,i,1,-1)=-bthet(2,i,1,1)
649 athet(1,i,-1,1)=-athet(1,i,1,1)
650 athet(2,i,-1,1)=-athet(2,i,1,1)
651 bthet(1,i,-1,1)=bthet(1,i,1,1)
652 bthet(2,i,-1,1)=bthet(2,i,1,1)
654 !elwrite(iout,*) "parmread kontrol in cryst_theta"
657 athet(1,i,-1,-1)=athet(1,-i,1,1)
658 athet(2,i,-1,-1)=-athet(2,-i,1,1)
659 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
660 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
661 athet(1,i,-1,1)=athet(1,-i,1,1)
662 athet(2,i,-1,1)=-athet(2,-i,1,1)
663 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
664 bthet(2,i,-1,1)=bthet(2,-i,1,1)
665 athet(1,i,1,-1)=-athet(1,-i,1,1)
666 athet(2,i,1,-1)=athet(2,-i,1,1)
667 bthet(1,i,1,-1)=bthet(1,-i,1,1)
668 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
673 polthet(j,i)=polthet(j,-i)
676 gthet(j,i)=gthet(j,-i)
679 !elwrite(iout,*) "parmread kontrol in cryst_theta"
681 !elwrite(iout,*) "parmread kontrol in cryst_theta"
684 ! & 'Parameters of the virtual-bond valence angles:'
685 ! write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
686 ! & ' ATHETA0 ',' A1 ',' A2 ',
689 ! write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
690 ! & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
692 ! write (iout,'(/a/9x,5a/79(1h-))')
693 ! & 'Parameters of the expression for sigma(theta_c):',
694 ! & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
695 ! & ' ALPH3 ',' SIGMA0C '
697 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
698 ! & (polthet(j,i),j=0,3),sigc0(i)
700 ! write (iout,'(/a/9x,5a/79(1h-))')
701 ! & 'Parameters of the second gaussian:',
702 ! & ' THETA0 ',' SIGMA0 ',' G1 ',
705 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
706 ! & sig0(i),(gthet(j,i),j=1,3)
709 'Parameters of the virtual-bond valence angles:'
710 write (iout,'(/a/9x,5a/79(1h-))') &
711 'Coefficients of expansion',&
712 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
713 ' b1*10^1 ',' b2*10^1 '
715 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
716 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
717 (10*bthet(j,i,1,1),j=1,2)
719 write (iout,'(/a/9x,5a/79(1h-))') &
720 'Parameters of the expression for sigma(theta_c):',&
721 ' alpha0 ',' alph1 ',' alph2 ',&
722 ' alhp3 ',' sigma0c '
724 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
725 (polthet(j,i),j=0,3),sigc0(i)
727 write (iout,'(/a/9x,5a/79(1h-))') &
728 'Parameters of the second gaussian:',&
729 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
732 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
733 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
738 ! Read the parameters of Utheta determined from ab initio surfaces
739 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
741 ! write (iout,*) "tu dochodze"
742 read (ithep,*) nthetyp,ntheterm,ntheterm2,&
743 ntheterm3,nsingle,ndouble
744 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
746 !----------------------------------------------------
747 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
748 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
749 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
750 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
751 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
752 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
753 !(maxtheterm,-maxthetyp1:maxthetyp1,&
754 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
755 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
756 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
757 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
758 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
759 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
760 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
761 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
762 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
763 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
764 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
765 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
766 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
767 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
768 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
769 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
770 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
773 read (ithep,*) (ithetyp(i),i=1,ntyp1)
775 ithetyp(i)=-ithetyp(-i)
777 ! write (iout,*) "tu dochodze"
778 aa0thet(:,:,:,:)=0.0d0
779 aathet(:,:,:,:,:)=0.0d0
780 bbthet(:,:,:,:,:,:)=0.0d0
781 ccthet(:,:,:,:,:,:)=0.0d0
782 ddthet(:,:,:,:,:,:)=0.0d0
783 eethet(:,:,:,:,:,:)=0.0d0
784 ffthet(:,:,:,:,:,:,:)=0.0d0
785 ggthet(:,:,:,:,:,:,:)=0.0d0
789 do j=-nthetyp,nthetyp
790 do k=-nthetyp,nthetyp
791 read (ithep,'(6a)') res1
792 read (ithep,*) aa0thet(i,j,k,iblock)
793 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
795 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
796 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
797 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
798 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
801 (((ffthet(llll,lll,ll,i,j,k,iblock),&
802 ffthet(lll,llll,ll,i,j,k,iblock),&
803 ggthet(llll,lll,ll,i,j,k,iblock),&
804 ggthet(lll,llll,ll,i,j,k,iblock),&
805 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
810 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
811 ! coefficients of theta-and-gamma-dependent terms are zero.
816 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
817 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
819 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
820 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
823 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
825 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
828 ! Substitution for D aminoacids from symmetry.
831 do j=-nthetyp,nthetyp
832 do k=-nthetyp,nthetyp
833 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
835 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
839 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
840 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
841 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
842 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
848 ffthet(llll,lll,ll,i,j,k,iblock)= &
849 ffthet(llll,lll,ll,-i,-j,-k,iblock)
850 ffthet(lll,llll,ll,i,j,k,iblock)= &
851 ffthet(lll,llll,ll,-i,-j,-k,iblock)
852 ggthet(llll,lll,ll,i,j,k,iblock)= &
853 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
854 ggthet(lll,llll,ll,i,j,k,iblock)= &
855 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
865 ! Control printout of the coefficients of virtual-bond-angle potentials
869 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
873 write (iout,'(//4a)') &
874 'Type ',onelett(i),onelett(j),onelett(k)
875 write (iout,'(//a,10x,a)') " l","a[l]"
876 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
877 write (iout,'(i2,1pe15.5)') &
878 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
880 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
881 "b",l,"c",l,"d",l,"e",l
883 write (iout,'(i2,4(1pe15.5))') m,&
884 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
885 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
889 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
890 "f+",l,"f-",l,"g+",l,"g-",l
893 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
894 ffthet(n,m,l,i,j,k,iblock),&
895 ffthet(m,n,l,i,j,k,iblock),&
896 ggthet(n,m,l,i,j,k,iblock),&
897 ggthet(m,n,l,i,j,k,iblock)
908 !-------------------------------------------
909 allocate(nlob(ntyp1)) !(ntyp1)
910 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
911 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
912 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
930 gaussc(l,k,j,i)=0.0D0
938 ! Read the parameters of the probability distribution/energy expression
939 ! of the side chains.
942 !c write (iout,*) "tu dochodze",i
943 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
947 dsc_inv(i)=1.0D0/dsc(i)
958 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
959 censc(1,1,-i)=censc(1,1,i)
960 censc(2,1,-i)=censc(2,1,i)
961 censc(3,1,-i)=-censc(3,1,i)
963 read (irotam,*) bsc(j,i)
964 read (irotam,*) (censc(k,j,i),k=1,3),&
965 ((blower(k,l,j),l=1,k),k=1,3)
966 censc(1,j,-i)=censc(1,j,i)
967 censc(2,j,-i)=censc(2,j,i)
968 censc(3,j,-i)=-censc(3,j,i)
969 ! BSC is amplitude of Gaussian
976 akl=akl+blower(k,m,j)*blower(l,m,j)
980 if (((k.eq.3).and.(l.ne.3)) &
981 .or.((l.eq.3).and.(k.ne.3))) then
982 gaussc(k,l,j,-i)=-akl
983 gaussc(l,k,j,-i)=-akl
995 write (iout,'(/a)') 'Parameters of side-chain local geometry'
999 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1000 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1001 ! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1002 ! write (iout,'(a,f10.4,4(16x,f10.4))')
1003 ! & 'Center ',(bsc(j,i),j=1,nlobi)
1004 ! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
1005 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1006 'log h',(bsc(j,i),j=1,nlobi)
1007 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1008 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1009 ! write (iout,'(a)')
1015 ! blower(k,l,j)=gaussc(ind,j,i)
1020 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1021 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1028 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1029 ! added by Urszula Kozlowska 07/11/2007
1031 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1039 read(irotam,*) sc_parmin(j,i)
1047 ! Read torsional parameters in old format
1049 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1051 read (itorp,*) ntortyp,nterm_old
1052 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1053 read (itorp,*) (itortyp(i),i=1,ntyp)
1055 !el from energy module--------
1056 allocate(v1(nterm_old,ntortyp,ntortyp))
1057 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1058 !el---------------------------
1064 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1070 write (iout,'(/a/)') 'Torsional constants:'
1073 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1074 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1082 ! Read torsional parameters
1084 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1086 read (itorp,*) ntortyp
1087 read (itorp,*) (itortyp(i),i=1,ntyp)
1088 write (iout,*) 'ntortyp',ntortyp
1090 !el from energy module---------
1091 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1092 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1094 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1095 allocate(vlor2(maxlor,ntortyp,ntortyp))
1096 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1097 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1099 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1100 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1101 !el---------------------------
1103 do i=-ntortyp,ntortyp
1104 do j=-ntortyp,ntortyp
1110 !el---------------------------
1114 itortyp(i)=-itortyp(-i)
1116 ! write (iout,*) 'ntortyp',ntortyp
1118 do j=-ntortyp+1,ntortyp-1
1119 read (itorp,*) nterm(i,j,iblock),&
1121 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1122 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1125 do k=1,nterm(i,j,iblock)
1126 read (itorp,*) kk,v1(k,i,j,iblock),&
1128 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1129 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1130 v0ij=v0ij+si*v1(k,i,j,iblock)
1133 do k=1,nlor(i,j,iblock)
1134 read (itorp,*) kk,vlor1(k,i,j),&
1135 vlor2(k,i,j),vlor3(k,i,j)
1136 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1139 v0(-i,-j,iblock)=v0ij
1146 write (iout,'(/a/)') 'Torsional constants:'
1149 write (iout,*) 'ityp',i,' jtyp',j
1150 write (iout,*) 'Fourier constants'
1151 do k=1,nterm(i,j,iblock)
1152 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1155 write (iout,*) 'Lorenz constants'
1156 do k=1,nlor(i,j,iblock)
1157 write (iout,'(3(1pe15.5))') &
1158 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1165 ! 6/23/01 Read parameters for double torsionals
1167 !el from energy module------------
1168 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1169 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1170 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1171 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1172 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1173 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1174 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1175 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1176 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1177 !---------------------------------
1181 do j=-ntortyp+1,ntortyp-1
1182 do k=-ntortyp+1,ntortyp-1
1183 read (itordp,'(3a1)') t1,t2,t3
1184 ! write (iout,*) "OK onelett",
1187 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1188 .or. t3.ne.toronelet(k)) then
1189 write (iout,*) "Error in double torsional parameter file",&
1192 call MPI_Finalize(Ierror)
1194 stop "Error in double torsional parameter file"
1196 read (itordp,*) ntermd_1(i,j,k,iblock),&
1197 ntermd_2(i,j,k,iblock)
1198 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1199 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1200 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1201 ntermd_1(i,j,k,iblock))
1202 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1203 ntermd_1(i,j,k,iblock))
1204 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1205 ntermd_1(i,j,k,iblock))
1206 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1207 ntermd_1(i,j,k,iblock))
1208 ! Martix of D parameters for one dimesional foureir series
1209 do l=1,ntermd_1(i,j,k,iblock)
1210 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1211 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1212 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1213 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1214 ! write(iout,*) "whcodze" ,
1215 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1217 read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1218 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1219 v2s(m,l,i,j,k,iblock),&
1220 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1221 ! Martix of D parameters for two dimesional fourier series
1222 do l=1,ntermd_2(i,j,k,iblock)
1224 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1225 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1226 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1227 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1236 write (iout,*) 'Constants for double torsionals'
1239 do j=-ntortyp+1,ntortyp-1
1240 do k=-ntortyp+1,ntortyp-1
1241 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1242 ' nsingle',ntermd_1(i,j,k,iblock),&
1243 ' ndouble',ntermd_2(i,j,k,iblock)
1245 write (iout,*) 'Single angles:'
1246 do l=1,ntermd_1(i,j,k,iblock)
1247 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1248 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1249 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1250 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1253 write (iout,*) 'Pairs of angles:'
1254 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1255 do l=1,ntermd_2(i,j,k,iblock)
1256 write (iout,'(i5,20f10.5)') &
1257 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1260 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1261 do l=1,ntermd_2(i,j,k,iblock)
1262 write (iout,'(i5,20f10.5)') &
1263 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1264 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1273 !elwrite(iout,*) "parmread kontrol sc-bb"
1274 ! Read of Side-chain backbone correlation parameters
1275 ! Modified 11 May 2012 by Adasko
1278 read (isccor,*) nsccortyp
1281 !c maxinter is maximum interaction sites
1282 !write(iout,*)"maxterm_sccor",maxterm_sccor
1283 !el from module energy-------------
1284 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1285 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1286 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1287 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1288 !-----------------------------------
1289 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1290 !-----------------------------------
1291 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1292 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1293 -nsccortyp:nsccortyp))
1294 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1295 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1296 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1297 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1298 !-----------------------------------
1299 do i=-nsccortyp,nsccortyp
1300 do j=-nsccortyp,nsccortyp
1304 !-----------------------------------
1306 read (isccor,*) (isccortyp(i),i=1,ntyp)
1308 isccortyp(i)=-isccortyp(-i)
1310 iscprol=isccortyp(20)
1311 ! write (iout,*) 'ntortyp',ntortyp
1313 !c maxinter is maximum interaction sites
1318 nterm_sccor(i,j),nlor_sccor(i,j)
1324 nterm_sccor(-i,j)=nterm_sccor(i,j)
1325 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1326 nterm_sccor(i,-j)=nterm_sccor(i,j)
1327 do k=1,nterm_sccor(i,j)
1328 read (isccor,*) kk,v1sccor(k,l,i,j),&
1330 if (j.eq.iscprol) then
1331 if (i.eq.isccortyp(10)) then
1332 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1333 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1335 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1336 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1337 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1338 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1339 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1340 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1341 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1342 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1345 if (i.eq.isccortyp(10)) then
1346 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1347 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1349 if (j.eq.isccortyp(10)) then
1350 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1351 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1353 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1354 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1355 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1356 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1357 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1358 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1362 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1363 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1364 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1365 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1368 do k=1,nlor_sccor(i,j)
1369 read (isccor,*) kk,vlor1sccor(k,i,j),&
1370 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1371 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1372 (1+vlor3sccor(k,i,j)**2)
1374 v0sccor(l,i,j)=v0ijsccor
1375 v0sccor(l,-i,j)=v0ijsccor1
1376 v0sccor(l,i,-j)=v0ijsccor2
1377 v0sccor(l,-i,-j)=v0ijsccor3
1383 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1386 write (iout,*) 'ityp',i,' jtyp',j
1387 write (iout,*) 'Fourier constants'
1388 do k=1,nterm_sccor(i,j)
1389 write (iout,'(2(1pe15.5))') &
1390 (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
1392 write (iout,*) 'Lorenz constants'
1393 do k=1,nlor_sccor(i,j)
1394 write (iout,'(3(1pe15.5))') &
1395 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1401 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1402 ! interaction energy of the Gly, Ala, and Pro prototypes.
1404 read (ifourier,*) nloctyp
1405 !el write(iout,*)"nloctyp",nloctyp
1406 !el from module energy-------
1407 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1408 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1409 allocate(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1410 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1411 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1412 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1413 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1414 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1416 do ii=-nloctyp-1,nloctyp+1
1424 ctilde(j,i,ii)=0.0d0
1425 dtilde(j,i,ii)=0.0d0
1429 !--------------------------------
1430 allocate(b(13,0:nloctyp))
1434 read (ifourier,*) (b(ii,i),ii=1,13)
1436 write (iout,*) 'Type',i
1437 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1445 B1tilde(1,i) = b(3,i)
1446 B1tilde(2,i) =-b(5,i)
1447 B1tilde(1,-i) =-b(3,i)
1448 B1tilde(2,-i) =b(5,i)
1449 ! b1tilde(1,i)=0.0d0
1450 ! b1tilde(2,i)=0.0d0
1470 Ctilde(1,1,i)=b(7,i)
1471 Ctilde(1,2,i)=b(9,i)
1472 Ctilde(2,1,i)=-b(9,i)
1473 Ctilde(2,2,i)=b(7,i)
1474 Ctilde(1,1,-i)=b(7,i)
1475 Ctilde(1,2,-i)=-b(9,i)
1476 Ctilde(2,1,-i)=b(9,i)
1477 Ctilde(2,2,-i)=b(7,i)
1479 ! Ctilde(1,1,i)=0.0d0
1480 ! Ctilde(1,2,i)=0.0d0
1481 ! Ctilde(2,1,i)=0.0d0
1482 ! Ctilde(2,2,i)=0.0d0
1495 Dtilde(1,1,i)=b(6,i)
1496 Dtilde(1,2,i)=b(8,i)
1497 Dtilde(2,1,i)=-b(8,i)
1498 Dtilde(2,2,i)=b(6,i)
1499 Dtilde(1,1,-i)=b(6,i)
1500 Dtilde(1,2,-i)=-b(8,i)
1501 Dtilde(2,1,-i)=b(8,i)
1502 Dtilde(2,2,-i)=b(6,i)
1504 ! Dtilde(1,1,i)=0.0d0
1505 ! Dtilde(1,2,i)=0.0d0
1506 ! Dtilde(2,1,i)=0.0d0
1507 ! Dtilde(2,2,i)=0.0d0
1508 EE(1,1,i)= b(10,i)+b(11,i)
1509 EE(2,2,i)=-b(10,i)+b(11,i)
1510 EE(2,1,i)= b(12,i)-b(13,i)
1511 EE(1,2,i)= b(12,i)+b(13,i)
1512 EE(1,1,-i)= b(10,i)+b(11,i)
1513 EE(2,2,-i)=-b(10,i)+b(11,i)
1514 EE(2,1,-i)=-b(12,i)+b(13,i)
1515 EE(1,2,-i)=-b(12,i)-b(13,i)
1521 ! ee(2,1,i)=ee(1,2,i)
1526 write (iout,*) 'Type',i
1528 ! write (iout,'(f10.5)') B1(:,i)
1529 write(iout,*) B1(1,i),B1(2,i)
1531 ! write (iout,'(f10.5)') B2(:,i)
1532 write(iout,*) B2(1,i),B2(2,i)
1535 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1539 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1543 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1548 ! Read electrostatic-interaction parameters
1551 write (iout,'(/a)') 'Electrostatic interaction constants:'
1552 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1553 'IT','JT','APP','BPP','AEL6','AEL3'
1555 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1556 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1557 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1558 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1563 app (i,j)=epp(i,j)*rri*rri
1564 bpp (i,j)=-2.0D0*epp(i,j)*rri
1565 ael6(i,j)=elpp6(i,j)*4.2D0**6
1566 ael3(i,j)=elpp3(i,j)*4.2D0**3
1567 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1572 ! Read side-chain interaction parameters.
1574 !el from module energy - COMMON.INTERACT-------
1575 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1576 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1577 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1578 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1579 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1580 allocate(epslip(ntyp,ntyp))
1591 !--------------------------------
1593 read (isidep,*) ipot,expon
1594 !el if (ipot.lt.1 .or. ipot.gt.5) then
1595 ! write (iout,'(2a)') 'Error while reading SC interaction',&
1596 ! 'potential file - unknown potential type.'
1600 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1601 ', exponents are ',expon,2*expon
1602 ! goto (10,20,30,30,40) ipot
1604 !----------------------- LJ potential ---------------------------------
1606 ! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1607 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1609 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1610 write (iout,'(a/)') 'The epsilon array:'
1611 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1612 write (iout,'(/a)') 'One-body parameters:'
1613 write (iout,'(a,4x,a)') 'residue','sigma'
1614 write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
1617 !----------------------- LJK potential --------------------------------
1619 ! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1620 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1621 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1623 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1624 write (iout,'(a/)') 'The epsilon array:'
1625 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1626 write (iout,'(/a)') 'One-body parameters:'
1627 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1628 write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
1632 !---------------------- GB or BP potential -----------------------------
1636 read (isidep,*)(eps(i,j),j=i,ntyp)
1638 read (isidep,*)(sigma0(i),i=1,ntyp)
1639 read (isidep,*)(sigii(i),i=1,ntyp)
1640 read (isidep,*)(chip(i),i=1,ntyp)
1641 read (isidep,*)(alp(i),i=1,ntyp)
1643 read (isidep,*)(epslip(i,j),j=i,ntyp)
1645 ! For the GB potential convert sigma'**2 into chi'
1648 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1652 write (iout,'(/a/)') 'Parameters of the BP potential:'
1653 write (iout,'(a/)') 'The epsilon array:'
1654 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1655 write (iout,'(/a)') 'One-body parameters:'
1656 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
1658 write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
1659 chip(i),alp(i),i=1,ntyp)
1662 !--------------------- GBV potential -----------------------------------
1664 ! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1665 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1666 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
1667 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1669 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1670 write (iout,'(a/)') 'The epsilon array:'
1671 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1672 write (iout,'(/a)') 'One-body parameters:'
1673 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
1674 's||/s_|_^2',' chip ',' alph '
1675 write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
1676 sigii(i),chip(i),alp(i),i=1,ntyp)
1679 write (iout,'(2a)') 'Error while reading SC interaction',&
1680 'potential file - unknown potential type.'
1686 !-----------------------------------------------------------------------
1687 ! Calculate the "working" parameters of SC interactions.
1689 !el from module energy - COMMON.INTERACT-------
1690 ! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
1691 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
1692 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
1693 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
1705 !--------------------------------
1710 epslip(i,j)=epslip(j,i)
1715 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1716 sigma(j,i)=sigma(i,j)
1717 rs0(i,j)=dwa16*sigma(i,j)
1721 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
1722 'Working parameters of the SC interactions:',&
1723 ' a ',' b ',' augm ',' sigma ',' r0 ',&
1728 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1737 sigeps=dsign(1.0D0,epsij)
1739 aa_aq(i,j)=epsij*rrij*rrij
1740 bb_aq(i,j)=-sigeps*epsij*rrij
1741 aa_aq(j,i)=aa_aq(i,j)
1742 bb_aq(j,i)=bb_aq(i,j)
1743 epsijlip=epslip(i,j)
1744 sigeps=dsign(1.0D0,epsijlip)
1745 epsijlip=dabs(epsijlip)
1746 aa_lip(i,j)=epsijlip*rrij*rrij
1747 bb_lip(i,j)=-sigeps*epsijlip*rrij
1748 aa_lip(j,i)=aa_lip(i,j)
1749 bb_lip(j,i)=bb_lip(i,j)
1751 sigt1sq=sigma0(i)**2
1752 sigt2sq=sigma0(j)**2
1755 ratsig1=sigt2sq/sigt1sq
1756 ratsig2=1.0D0/ratsig1
1757 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1758 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1759 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1763 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1764 sigmaii(i,j)=rsum_max
1765 sigmaii(j,i)=rsum_max
1767 ! sigmaii(i,j)=r0(i,j)
1768 ! sigmaii(j,i)=r0(i,j)
1770 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1771 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1772 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1773 augm(i,j)=epsij*r_augm**(2*expon)
1774 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1781 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
1782 restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
1783 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1788 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
1796 ! Define the SC-p interaction constants
1806 !elwrite(iout,*) "parmread kontrol before oldscp"
1808 ! Define the SC-p interaction constants
1812 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1814 ! aad(i,1)=0.3D0*4.0D0**12
1815 ! Following line for constants currently implemented
1816 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
1817 aad(i,1)=1.5D0*4.0D0**12
1818 ! aad(i,1)=0.17D0*5.6D0**12
1820 ! "Soft" SC-p repulsion
1822 ! Following line for constants currently implemented
1823 ! aad(i,1)=0.3D0*4.0D0**6
1824 ! "Hard" SC-p repulsion
1825 bad(i,1)=3.0D0*4.0D0**6
1826 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1835 ! 8/9/01 Read the SC-p interaction constants from file
1838 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1841 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1842 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1843 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1844 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1848 write (iout,*) "Parameters of SC-p interactions:"
1850 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
1851 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1856 ! Define the constants of the disulfide bridge
1860 ! Old arbitrary potential - commented out.
1865 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1866 ! energy surface of diethyl disulfide.
1867 ! A. Liwo and U. Kozlowska, 11/24/03
1878 write (iout,'(/a)') "Disulfide bridge parameters:"
1879 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1880 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1881 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1882 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
1886 end subroutine parmread
1888 !-----------------------------------------------------------------------------
1890 !-----------------------------------------------------------------------------
1891 subroutine mygetenv(string,var)
1895 ! This subroutine passes the environmental variables to FORTRAN program.
1896 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
1897 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
1898 ! reads the environmental variables from $HOME/.env
1900 ! Usage: As for the standard FORTRAN GETENV subroutine.
1902 ! Purpose: some versions/installations of MPI do not transfer the environmental
1903 ! variables to slave processors, if these variables are set in the shell script
1904 ! from which mpirun is called.
1913 character*(*) :: string,var
1914 #if defined(MYGETENV) && defined(MPI)
1915 ! include "DIMENSIONS.ZSCOPT"
1917 ! include "COMMON.MPI"
1918 !el character*360 ucase
1920 character(len=360) :: string1(360),karta
1921 character(len=240) :: home
1924 call getenv("HOME",home)
1925 open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
1927 read (99,end=111,err=111,'(a)') karta
1931 call split_string(karta,string1,80,n)
1932 if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
1933 string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
1935 print *,"Processor",me,": ",var(:ilen(var)),&
1936 " assigned to ",string(:ilen(string))
1941 111 print *,"Environment variable ",string(:ilen(string))," not set."
1944 112 print *,"Error opening environment file!"
1946 call getenv(string,var)
1949 end subroutine mygetenv
1950 !-----------------------------------------------------------------------------
1952 !-----------------------------------------------------------------------------
1953 subroutine read_general_data(*)
1955 use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions
1956 use energy_data, only:distchainmax
1957 use geometry_data, only:boxxsize,boxysize,boxzsize
1959 ! include "DIMENSIONS"
1960 ! include "DIMENSIONS.ZSCOPT"
1961 ! include "DIMENSIONS.FREE"
1962 ! include "COMMON.TORSION"
1963 ! include "COMMON.INTERACT"
1964 ! include "COMMON.IOUNITS"
1965 ! include "COMMON.TIME1"
1966 ! include "COMMON.PROT"
1967 ! include "COMMON.PROTFILES"
1968 ! include "COMMON.CHAIN"
1969 ! include "COMMON.NAMES"
1970 ! include "COMMON.FFIELD"
1971 ! include "COMMON.ENEPS"
1972 ! include "COMMON.WEIGHTS"
1973 ! include "COMMON.FREE"
1974 ! include "COMMON.CONTROL"
1975 ! include "COMMON.ENERGIES"
1976 character(len=800) :: controlcard
1977 integer :: i,j,k,ii,n_ene_found
1978 integer :: ind,itype1,itype2,itypf,itypsc,itypp
1981 !el character*16 ucase
1982 character(len=16) :: key
1984 call card_concat(controlcard,.true.)
1985 call readi(controlcard,"N_ENE",n_eneW,max_eneW)
1986 if (n_eneW.gt.max_eneW) then
1987 write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
1991 call readi(controlcard,"NPARMSET",nparmset,1)
1992 !elwrite(iout,*)"in read_gen data"
1993 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
1994 call readi(controlcard,"IPARMPRINT",iparmprint,1)
1995 write (iout,*) "PARMPRINT",iparmprint
1996 if (nparmset.gt.max_parm) then
1997 write (iout,*) "Error: parameter out of range: NPARMSET",&
2001 !elwrite(iout,*)"in read_gen data"
2002 call readi(controlcard,"MAXIT",maxit,5000)
2003 call reada(controlcard,"FIMIN",fimin,1.0d-3)
2004 call readi(controlcard,"ENSEMBLES",ensembles,0)
2005 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
2006 write (iout,*) "Number of energy parameter sets",nparmset
2007 allocate(isampl(nparmset))
2008 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
2009 write (iout,*) "MaxSlice",MaxSlice
2010 call readi(controlcard,"NSLICE",nslice,1)
2011 !elwrite(iout,*)"in read_gen data"
2013 if (nslice.gt.MaxSlice) then
2014 write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
2018 write (iout,*) "Frequency of storing conformations",&
2019 (isampl(i),i=1,nparmset)
2020 write (iout,*) "Maxit",maxit," Fimin",fimin
2021 call readi(controlcard,"NQ",nQ,1)
2022 if (nQ.gt.MaxQ) then
2023 write (iout,*) "Error: parameter out of range: NQ",nq,&
2028 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
2029 call reada(controlcard,"DELTA",delta,1.0d-2)
2030 call readi(controlcard,"EINICHECK",einicheck,2)
2031 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
2032 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
2033 call readi(controlcard,"RESCALE",rescale_modeW,1)
2034 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2035 call reada(controlcard,'BOXY',boxysize,100.0d0)
2036 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2037 ions=index(controlcard,"IONS").gt.0
2038 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2039 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2040 write(iout,*) "R_CUT_ELE=",r_cut_ele
2041 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
2042 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
2043 call readi(controlcard,'SYM',symetr,1)
2044 write (iout,*) "DISTCHAINMAX",distchainmax
2045 write (iout,*) "delta",delta
2046 write (iout,*) "einicheck",einicheck
2047 write (iout,*) "rescale_mode",rescale_modeW
2049 bxfile=index(controlcard,"BXFILE").gt.0
2050 cxfile=index(controlcard,"CXFILE").gt.0
2051 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
2053 histfile=index(controlcard,"HISTFILE").gt.0
2054 histout=index(controlcard,"HISTOUT").gt.0
2055 entfile=index(controlcard,"ENTFILE").gt.0
2056 zscfile=index(controlcard,"ZSCFILE").gt.0
2057 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
2058 write (iout,*) "with_dihed_constr ",with_dihed_constr
2059 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2061 end subroutine read_general_data
2062 !------------------------------------------------------------------------------
2063 subroutine read_efree(*)
2065 ! Read molecular data
2068 ! include 'DIMENSIONS'
2069 ! include 'DIMENSIONS.ZSCOPT'
2070 ! include 'DIMENSIONS.COMPAR'
2071 ! include 'DIMENSIONS.FREE'
2072 ! include 'COMMON.IOUNITS'
2073 ! include 'COMMON.TIME1'
2074 ! include 'COMMON.SBRIDGE'
2075 ! include 'COMMON.CONTROL'
2076 ! include 'COMMON.CHAIN'
2077 ! include 'COMMON.HEADER'
2078 ! include 'COMMON.GEO'
2079 ! include 'COMMON.FREE'
2080 character(len=320) :: controlcard !,ucase
2081 integer :: iparm,ib,i,j,npars
2091 ! call alloc_wham_arrays
2092 ! allocate(nT_h(nParmSet))
2093 ! allocate(replica(nParmSet))
2094 ! allocate(umbrella(nParmSet))
2095 ! allocate(read_iset(nParmSet))
2096 ! allocate(nT_h(nParmSet))
2100 call card_concat(controlcard,.true.)
2101 call readi(controlcard,'NT',nT_h(iparm),1)
2102 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
2104 if (nT_h(iparm).gt.MaxT_h) then
2105 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
2109 replica(iparm)=index(controlcard,"REPLICA").gt.0
2110 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
2111 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
2112 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
2113 replica(iparm)," umbrella ",umbrella(iparm),&
2114 " read_iset",read_iset(iparm)
2117 call card_concat(controlcard,.true.)
2118 call readi(controlcard,'NR',nR(ib,iparm),1)
2119 if (umbrella(iparm)) then
2122 nRR(ib,iparm)=nR(ib,iparm)
2124 if (nR(ib,iparm).gt.MaxR) then
2125 write (iout,*) "Error: parameter out of range: NR",&
2129 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
2130 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
2131 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
2134 call card_concat(controlcard,.true.)
2135 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
2137 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
2142 write (iout,*) "ib",ib," beta_h",&
2143 1.0d0/(0.001987*beta_h(ib,iparm))
2144 write (iout,*) "nR",nR(ib,iparm)
2145 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
2147 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
2148 "q0",(q0(j,i,ib,iparm),j=1,nQ)
2160 nR(ib,iparm)=nR(ib,1)
2161 if (umbrella(iparm)) then
2164 nRR(ib,iparm)=nR(ib,1)
2166 beta_h(ib,iparm)=beta_h(ib,1)
2168 f(i,ib,iparm)=f(i,ib,1)
2170 KH(j,i,ib,iparm)=KH(j,i,ib,1)
2171 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
2174 replica(iparm)=replica(1)
2175 umbrella(iparm)=umbrella(1)
2176 read_iset(iparm)=read_iset(1)
2183 end subroutine read_efree
2184 !-----------------------------------------------------------------------------
2185 subroutine read_protein_data(*)
2187 ! include "DIMENSIONS"
2188 ! include "DIMENSIONS.ZSCOPT"
2189 ! include "DIMENSIONS.FREE"
2193 integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
2194 ! include "COMMON.MPI"
2196 ! include "COMMON.CHAIN"
2197 ! include "COMMON.IOUNITS"
2198 ! include "COMMON.PROT"
2199 ! include "COMMON.PROTFILES"
2200 ! include "COMMON.NAMES"
2201 ! include "COMMON.FREE"
2202 ! include "COMMON.OBCINKA"
2203 character(len=64) :: nazwa
2204 character(len=16000) :: controlcard
2205 integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
2206 !el external ilen,iroof
2215 ! Read names of files with conformation data.
2216 if (replica(iparm)) then
2222 do ii=1,nRR(ib,iparm)
2223 write (iout,*) "Parameter set",iparm," temperature",ib,&
2226 call card_concat(controlcard,.true.)
2227 write (iout,*) controlcard(:ilen(controlcard))
2228 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
2229 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
2230 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
2231 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
2232 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
2233 maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
2234 call reada(controlcard,"TIME_START",&
2235 time_start_collect(ii,ib,iparm),0.0d0)
2236 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
2238 write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
2239 " rec_end",rec_end(ii,ib,iparm)
2240 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
2241 " time_end",time_end_collect(ii,ib,iparm)
2243 if (replica(iparm)) then
2244 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
2245 write (iout,*) "Number of trajectories",totraj(ii,iparm)
2248 if (nfile_bin(ii,ib,iparm).lt.2 &
2249 .and. nfile_asc(ii,ib,iparm).eq.0 &
2250 .and. nfile_cx(ii,ib,iparm).eq.0) then
2251 write (iout,*) "Error - no action specified!"
2254 if (nfile_bin(ii,ib,iparm).gt.0) then
2255 call card_concat(controlcard,.false.)
2256 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
2257 maxfile_prot,nfile_bin(ii,ib,iparm))
2259 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
2260 write(iout,*) (protfiles(i,1,ii,ib,iparm),&
2261 i=1,nfile_bin(ii,ib,iparm))
2264 if (nfile_asc(ii,ib,iparm).gt.0) then
2265 call card_concat(controlcard,.false.)
2266 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2267 maxfile_prot,nfile_asc(ii,ib,iparm))
2269 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
2270 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2271 i=1,nfile_asc(ii,ib,iparm))
2273 else if (nfile_cx(ii,ib,iparm).gt.0) then
2274 call card_concat(controlcard,.false.)
2275 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2276 maxfile_prot,nfile_cx(ii,ib,iparm))
2278 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
2279 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2280 i=1,nfile_cx(ii,ib,iparm))
2289 end subroutine read_protein_data
2290 !-------------------------------------------------------------------------------
2291 subroutine readsss(rekord,lancuch,wartosc,default)
2293 character*(*) :: rekord,lancuch,wartosc,default
2294 character(len=80) :: aux
2295 integer :: lenlan,lenrec,iread,ireade
2299 lenlan=ilen(lancuch)
2301 iread=index(rekord,lancuch(:lenlan)//"=")
2302 ! print *,"rekord",rekord," lancuch",lancuch
2303 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2304 if (iread.eq.0) then
2308 iread=iread+lenlan+1
2309 ! print *,"iread",iread
2310 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2311 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2313 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2315 ! print *,"iread",iread
2316 if (iread.gt.lenrec) then
2321 ! print *,"ireade",ireade
2322 do while (ireade.lt.lenrec .and. &
2323 .not.iblnk(rekord(ireade:ireade)))
2326 wartosc=rekord(iread:ireade)
2328 end subroutine readsss
2329 !----------------------------------------------------------------------------
2330 subroutine multreads(rekord,lancuch,tablica,dim,default)
2333 character*(*) rekord,lancuch,tablica(dim),default
2334 character(len=80) :: aux
2335 integer :: lenlan,lenrec,iread,ireade
2342 lenlan=ilen(lancuch)
2344 iread=index(rekord,lancuch(:lenlan)//"=")
2345 ! print *,"rekord",rekord," lancuch",lancuch
2346 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2347 if (iread.eq.0) return
2348 iread=iread+lenlan+1
2350 ! print *,"iread",iread
2351 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2352 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2354 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2356 ! print *,"iread",iread
2357 if (iread.gt.lenrec) return
2359 ! print *,"ireade",ireade
2360 do while (ireade.lt.lenrec .and. &
2361 .not.iblnk(rekord(ireade:ireade)))
2364 tablica(i)=rekord(iread:ireade)
2367 end subroutine multreads
2368 !----------------------------------------------------------------------------
2369 subroutine split_string(rekord,tablica,dim,nsub)
2371 integer :: dim,nsub,i,ii,ll,kk
2372 character*(*) tablica(dim)
2373 character*(*) rekord
2383 ! Find the start of term name
2385 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2388 ! Parse the name into TABLICA(i) until blank found
2389 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2391 tablica(i)(kk:kk)=rekord(ii:ii)
2394 if (kk.gt.0) nsub=nsub+1
2395 if (ii.gt.ll) return
2398 end subroutine split_string
2399 !--------------------------------------------------------------------------------
2401 !--------------------------------------------------------------------------------
2402 subroutine read_compar
2404 ! Read molecular data
2406 use conform_compar, only:alloc_compar_arrays
2407 use control_data, only:pdbref
2408 use geometry_data, only:deg2rad,rad2deg
2410 ! include 'DIMENSIONS'
2411 ! include 'DIMENSIONS.ZSCOPT'
2412 ! include 'DIMENSIONS.COMPAR'
2413 ! include 'DIMENSIONS.FREE'
2414 ! include 'COMMON.IOUNITS'
2415 ! include 'COMMON.TIME1'
2416 ! include 'COMMON.SBRIDGE'
2417 ! include 'COMMON.CONTROL'
2418 ! include 'COMMON.COMPAR'
2419 ! include 'COMMON.CHAIN'
2420 ! include 'COMMON.HEADER'
2421 ! include 'COMMON.GEO'
2422 ! include 'COMMON.FREE'
2423 character(len=320) :: controlcard !,ucase
2424 character(len=64) :: wfile
2428 !elwrite(iout,*)"jestesmy w read_compar"
2429 call card_concat(controlcard,.true.)
2430 pdbref=(index(controlcard,'PDBREF').gt.0)
2431 call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
2432 call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
2433 call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
2434 call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
2435 verbose = index(controlcard,"VERBOSE").gt.0
2436 lgrp=index(controlcard,"STATIN").gt.0
2437 lgrp_out=index(controlcard,"STATOUT").gt.0
2438 merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
2439 binary = index(controlcard,"BINARY").gt.0
2440 rmscut_base_up=rmscut_base_up/50
2441 rmscut_base_low=rmscut_base_low/50
2442 call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
2443 call readi(controlcard,'NLEVEL',nlevel,1)
2444 if (nlevel.lt.0) then
2446 call alloc_compar_arrays(maxfrag,1)
2449 allocate(nfrag(nlevel))
2451 ! Read the data pertaining to elementary fragments (level 1)
2452 call readi(controlcard,'NFRAG',nfrag(1),0)
2453 write(iout,*)"nfrag(1)",nfrag(1)
2454 call alloc_compar_arrays(nfrag(1),nlevel)
2456 call card_concat(controlcard,.true.)
2457 write (iout,*) controlcard(:ilen(controlcard))
2458 call readi(controlcard,'NPIECE',npiece(j,1),0)
2459 call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
2460 call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
2461 call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
2462 call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
2463 call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
2464 call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
2465 call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
2466 call readi(controlcard,'RMS',irms(j,1),0)
2467 call readi(controlcard,'LOCAL',iloc(j),1)
2468 call readi(controlcard,'ELCONT',ielecont(j,1),1)
2469 if (ielecont(j,1).eq.0) then
2470 call readi(controlcard,'SCCONT',isccont(j,1),1)
2472 ang_cut(j)=ang_cut(j)*deg2rad
2473 ang_cut1(j)=ang_cut1(j)*deg2rad
2475 call card_concat(controlcard,.true.)
2476 call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
2477 call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
2479 write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
2480 (ifrag(1,k,j),ifrag(2,k,j),&
2481 k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
2482 " ang_cut1",ang_cut1(j)*rad2deg
2483 write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
2484 write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
2485 write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
2486 " ilocal",iloc(j)," isccont",isccont(j,1)
2488 ! Read data pertaning to higher levels
2490 call card_concat(controlcard,.true.)
2491 call readi(controlcard,'NFRAG',NFRAG(i),0)
2492 write (iout,*) "i",i," nfrag",nfrag(i)
2494 call card_concat(controlcard,.true.)
2496 call readi(controlcard,'ELCONT',ielecont(j,i),0)
2497 if (ielecont(j,i).eq.0) then
2498 call readi(controlcard,'SCCONT',isccont(j,i),1)
2500 call readi(controlcard,'RMS',irms(j,i),0)
2506 call readi(controlcard,'NPIECE',npiece(j,i),0)
2507 call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
2508 call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
2509 call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
2511 call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
2512 call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
2513 write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
2514 n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
2515 " isccont",isccont(j,i)," irms",irms(j,i)
2516 write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
2517 write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
2518 write(iout,*)"nc_frac",nc_fragm(j,i),&
2519 " nc_req",nc_req_setf(j,i)
2522 if (binary) write (iout,*) "Classes written in binary format."
2525 call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
2526 call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
2527 call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
2528 call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
2529 call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
2530 call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
2531 call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
2532 call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
2533 call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
2534 call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
2535 call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
2536 call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
2537 call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
2538 call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
2539 call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
2540 call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
2541 call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
2542 call readi(controlcard,'RMS_SINGLE',irms_single,0)
2543 call readi(controlcard,'CONT_SINGLE',icont_single,1)
2544 call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
2545 call readi(controlcard,'RMS_PAIR',irms_pair,0)
2546 call readi(controlcard,'CONT_PAIR',icont_pair,1)
2547 call readi(controlcard,'SPLIT_BET',isplit_bet,0)
2548 angcut_hel=angcut_hel*deg2rad
2549 angcut1_hel=angcut1_hel*deg2rad
2550 angcut_bet=angcut_bet*deg2rad
2551 angcut1_bet=angcut1_bet*deg2rad
2552 angcut_strand=angcut_strand*deg2rad
2553 angcut1_strand=angcut1_strand*deg2rad
2554 write (iout,*) "Automatic detection of structural elements"
2555 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2556 ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2557 ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
2558 ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2559 ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
2560 ' SPLIT_BET',isplit_bet
2561 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2562 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2563 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2564 ' MAXANG_HEL',angcut1_hel*rad2deg
2565 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2566 ' MAXANG_BET',angcut1_bet*rad2deg
2567 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2568 ' MAXANG_STRAND',angcut1_strand*rad2deg
2569 write (iout,*) 'FRAC_MIN',frac_min_set
2571 end subroutine read_compar
2572 !--------------------------------------------------------------------------------
2574 !--------------------------------------------------------------------------------
2575 subroutine read_ref_structure(*)
2577 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
2580 use control_data, only:pdbref
2581 use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
2582 nstart_sup,nstart_seq,nperm,nres0
2583 use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
2584 use compare, only:seq_comp !,contact,elecont
2585 use geometry, only:chainbuild,dist
2586 use io_config, only:readpdb
2588 use conform_compar, only:contact,elecont
2590 ! include 'DIMENSIONS'
2591 ! include 'DIMENSIONS.ZSCOPT'
2592 ! include 'DIMENSIONS.COMPAR'
2593 ! include 'COMMON.IOUNITS'
2594 ! include 'COMMON.GEO'
2595 ! include 'COMMON.VAR'
2596 ! include 'COMMON.INTERACT'
2597 ! include 'COMMON.LOCAL'
2598 ! include 'COMMON.NAMES'
2599 ! include 'COMMON.CHAIN'
2600 ! include 'COMMON.FFIELD'
2601 ! include 'COMMON.SBRIDGE'
2602 ! include 'COMMON.HEADER'
2603 ! include 'COMMON.CONTROL'
2604 ! include 'COMMON.CONTACTS1'
2605 ! include 'COMMON.PEPTCONT'
2606 ! include 'COMMON.TIME1'
2607 ! include 'COMMON.COMPAR'
2608 character(len=4) :: sequence(nres)
2610 !el real(kind=8) :: x(maxvar)
2611 integer :: itype_pdb(nres,5)
2612 !el logical seq_comp
2613 integer :: i,j,k,nres_pdb,iaux,mnum
2614 real(kind=8) :: ddsc !el,dist
2615 integer :: kkk !,ilen
2619 write (iout,*) "pdbref",pdbref
2621 read(inp,'(a)') pdbfile
2622 write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
2623 pdbfile(:ilen(pdbfile))
2624 open(ipdbin,file=pdbfile,status='old',err=33)
2626 33 write (iout,'(a)') 'Error opening PDB file.'
2631 itype_pdb(i,mnum)=itype(i,mnum)
2637 iaux=itype_pdb(i,mnum)
2638 itype_pdb(i,mnum)=itype(i,mnum)
2646 if (nsup.le.(nct-nnt+1)) then
2647 do i=0,nct-nnt+1-nsup
2648 if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
2650 do j=nnt+nsup-1,nnt,-1
2652 cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
2655 do j=nnt+nsup-1,nnt,-1
2657 cref(k,j+i,kkk)=cref(k,j,kkk)
2659 phi_ref(j+i)=phi_ref(j)
2660 theta_ref(j+i)=theta_ref(j)
2661 alph_ref(j+i)=alph_ref(j)
2662 omeg_ref(j+i)=omeg_ref(j)
2666 write (iout,'(i5,3f10.5,5x,3f10.5)') &
2667 j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
2675 write (iout,'(a)') &
2676 'Error - sequences to be superposed do not match.'
2679 do i=0,nsup-(nct-nnt+1)
2680 if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
2683 nstart_sup=nstart_sup+i
2688 write (iout,'(a)') &
2689 'Error - sequences to be superposed do not match.'
2693 write (iout,'(a,i5)') &
2694 'Experimental structure begins at residue',nstart_seq
2696 call read_angles(inp,*38)
2698 38 write (iout,'(a)') 'Error reading reference structure.'
2707 cref(j,i,kkk)=c(j,i)
2711 nend_sup=nstart_sup+nsup-1
2714 c(j,i)=cref(j,i,kkk)
2720 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
2722 if (itype(i,mnum).ne.10) then
2723 ddsc = dist(i,nres+i)
2725 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
2729 dc_norm(j,nres+i)=0.0d0
2732 ! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
2733 ! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
2734 ! dc_norm(3,nres+i)**2
2736 dc(j,i)=c(j,i+1)-c(j,i)
2740 dc_norm(j,i)=dc(j,i)/ddsc
2743 ! print *,"Calling contact"
2744 call contact(.true.,ncont_ref,icont_ref(1,1),&
2745 nstart_sup,nend_sup)
2746 ! print *,"Calling elecont"
2747 call elecont(.true.,ncont_pept_ref,&
2748 icont_pept_ref(1,1),&
2749 nstart_sup,nend_sup)
2750 write (iout,'(a,i3,a,i3,a,i3,a)') &
2751 'Number of residues to be superposed:',nsup,&
2752 ' (from residue',nstart_sup,' to residue',&
2755 end subroutine read_ref_structure
2756 !--------------------------------------------------------------------------------
2758 !--------------------------------------------------------------------------------
2759 subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
2761 use geometry_data, only:nres,c
2762 use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
2763 ! implicit real*8 (a-h,o-z)
2764 ! include 'DIMENSIONS'
2765 ! include 'DIMENSIONS.ZSCOPT'
2766 ! include 'COMMON.CHAIN'
2767 ! include 'COMMON.INTERACT'
2768 ! include 'COMMON.NAMES'
2769 ! include 'COMMON.IOUNITS'
2770 ! include 'COMMON.HEADER'
2771 ! include 'COMMON.SBRIDGE'
2772 character(len=50) :: tytul
2773 character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
2774 'D','E','F','G','H','I','J'/),shape(chainid))
2775 integer,dimension(nres) :: ica !(maxres)
2776 real(kind=8) :: temp,efree,etot,entropy,rmsdev
2777 integer :: ii,i,j,iti,ires,iatom,ichain,mnum
2778 write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
2780 write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
2782 write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
2790 if (iti.eq.ntyp1) then
2793 write (ipdb,'(a)') 'TER'
2798 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
2802 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
2803 ires,(c(j,nres+i),j=1,3)
2807 write (ipdb,'(a)') 'TER'
2810 if (itype(i,mnum).eq.ntyp1) cycle
2811 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
2812 write (ipdb,30) ica(i),ica(i+1)
2813 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
2814 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2815 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
2816 write (ipdb,30) ica(i),ica(i)+1
2819 if (itype(nct,molnum(nct)).ne.10) then
2820 write (ipdb,30) ica(nct),ica(nct)+1
2823 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2825 write (ipdb,'(a6)') 'ENDMDL'
2826 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2827 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2828 30 FORMAT ('CONECT',8I5)
2830 end subroutine pdboutW
2832 !------------------------------------------------------------------------------
2834 !-----------------------------------------------------------------------------
2835 !-----------------------------------------------------------------------------