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 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,5) :: itype_pdb !(maxres)
158 integer :: i,j,i1,i2,it1,it2,mnum
159 real(kind=8) :: scalscp
160 !el logical :: seq_comp
161 write(iout,*) "START MOLREAD"
165 print *,"nres_molec, initial",nres_molec(i)
168 call card_concat(controlcard,.true.)
169 call reada(controlcard,'SCAL14',scal14,0.4d0)
170 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
171 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
172 call reada(controlcard,'TEMP0',temp0,300.0d0) !el
173 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
174 r0_corr=cutoff_corr-delt_corr
175 call readi(controlcard,"NRES",nres_molec(1),0)
176 call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
177 call readi(controlcard,"NRES_CAT",nres_molec(5),0)
180 nres=nres_molec(i)+nres
182 ! allocate(sequence(nres+1))
183 !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
184 call alloc_geo_arrays
185 call alloc_ener_arrays
186 ! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
187 !----------------------------
188 allocate(c(3,2*nres+2))
189 allocate(dc(3,0:2*nres+2))
190 allocate(itype(nres+2,5))
191 allocate(itel(nres+2))
192 if (.not. allocated(molnum)) allocate(molnum(nres+2))
208 !--------------------------
210 allocate(sequence(nres+1))
211 iscode=index(controlcard,"ONE_LETTER")
213 write (iout,*) "Error: no residues in molecule"
216 if (nres.gt.maxres) then
217 write (iout,*) "Error: too many residues",nres,maxres
219 write(iout,*) 'nres=',nres
221 ! Read sequence of the protein
222 if (iscode.gt.0) then
223 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
225 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
227 ! Convert sequence to numeric code
231 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
233 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
236 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
238 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
241 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
244 write (iout,*) "Numeric code:"
245 write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres)
249 if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
251 if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
255 else if (iabs(itype(i+1,mnum)).ne.20) then
257 else if (iabs(itype(i,mnum)).ne.20) then
264 write (iout,*) "ITEL"
267 write (iout,*) i,itype(i,mnum),itel(i)
272 if (with_dihed_constr) then
274 read (inp,*) ndih_constr
275 if (ndih_constr.gt.0) then
277 write (iout,*) 'FTORS',ftors
278 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
280 'There are',ndih_constr,' constraints on phi angles.'
282 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
285 phi0(i)=deg2rad*phi0(i)
286 drange(i)=deg2rad*drange(i)
294 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
295 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
296 write(iout,*) 'NNT=',NNT,' NCT=',NCT
300 write (iout,'(/a,i3,a)') 'The chain contains',ns,&
301 ' disulfide-bridging cysteines.'
302 write (iout,'(20i4)') (iss(i),i=1,ns)
303 write (iout,'(/a/)') 'Pre-formed links are:'
310 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
311 restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',&
312 dhpb(i),ebr,forcon(i)
317 end subroutine molread
318 !-----------------------------------------------------------------------------
320 !-----------------------------------------------------------------------------
321 subroutine parmread(iparm,*)
326 ! Read the parameters of the probability distributions of the virtual-bond
327 ! valence angles and the side chains and energy parameters.
333 use control_data, only: maxterm,maxlor,maxterm_sccor,& !maxtor
334 maxtermd_1,maxtermd_2 !,maxthetyp,maxthetyp1
338 use io_config, only: printmat
339 use control, only: getenv_loc
346 ! implicit real*8 (a-h,o-z)
347 ! include 'DIMENSIONS'
348 ! include 'DIMENSIONS.ZSCOPT'
349 ! include 'DIMENSIONS.FREE'
350 ! include 'COMMON.IOUNITS'
351 ! include 'COMMON.CHAIN'
352 ! include 'COMMON.INTERACT'
353 ! include 'COMMON.GEO'
354 ! include 'COMMON.LOCAL'
355 ! include 'COMMON.TORSION'
356 ! include 'COMMON.FFIELD'
357 ! include 'COMMON.NAMES'
358 ! include 'COMMON.SBRIDGE'
359 ! include 'COMMON.WEIGHTS'
360 ! include 'COMMON.ENEPS'
361 ! include 'COMMON.SCCOR'
362 ! include 'COMMON.SCROT'
363 ! include 'COMMON.FREE'
364 character(len=1) :: t1,t2,t3
365 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
366 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
368 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
369 character(len=800) :: controlcard
370 character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
371 tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
375 character(len=16) :: key
377 !el real(kind=8) :: ip,mp
378 real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
379 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
380 real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
382 integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n
383 integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
387 ! Set LPRINT=.TRUE. for debugging
388 dwa16=2.0d0**(1.0d0/6.0d0)
391 ! Assign virtual-bond length
394 vblinv2=vblinv*vblinv
396 call card_concat(controlcard,.true.)
399 allocate(ww(max_eneW))
401 key = wname(i)(:ilen(wname(i)))
402 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
405 write (iout,*) "iparm",iparm," myparm",myparm
406 ! If reading not own parameters, skip assignment
408 if (iparm.eq.myparm .or. .not.separate_parset) then
411 ! Setup weights for UNRES
436 allocate(weights(n_ene))
451 weights(15)=0 !wstrain !
452 weights(16)=0 !wvdwpp !
454 weights(18)=0 !scal14 !
460 call card_concat(controlcard,.false.)
462 ! Return if not own parameters
464 if (iparm.ne.myparm .and. separate_parset) return
466 call reads(controlcard,"BONDPAR",bondname_t,bondname)
467 open (ibond,file=bondname_t,status='old')
469 call reads(controlcard,"THETPAR",thetname_t,thetname)
470 open (ithep,file=thetname_t,status='old')
472 call reads(controlcard,"ROTPAR",rotname_t,rotname)
473 open (irotam,file=rotname_t,status='old')
475 call reads(controlcard,"TORPAR",torname_t,torname)
476 open (itorp,file=torname_t,status='old')
478 call reads(controlcard,"TORDPAR",tordname_t,tordname)
479 open (itordp,file=tordname_t,status='old')
481 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
482 open (isccor,file=sccorname_t,status='old')
484 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
485 open (ifourier,file=fouriername_t,status='old')
487 call reads(controlcard,"ELEPAR",elename_t,elename)
488 open (ielep,file=elename_t,status='old')
490 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
491 open (isidep,file=sidename_t,status='old')
493 call reads(controlcard,"SCPPAR",scpname_t,scpname)
494 open (iscpp,file=scpname_t,status='old')
496 call getenv_loc('IONPAR',ionname)
497 open (iion,file=ionname,status='old')
499 write (iout,*) "Parameter set:",iparm
500 write (iout,*) "Energy-term weights:"
502 write (iout,'(a16,f10.5)') wname(i),ww(i)
504 write (iout,*) "Sidechain potential file : ",&
505 sidename_t(:ilen(sidename_t))
507 write (iout,*) "SCp potential file : ",&
508 scpname_t(:ilen(scpname_t))
510 write (iout,*) "Electrostatic potential file : ",&
511 elename_t(:ilen(elename_t))
512 write (iout,*) "Cumulant coefficient file : ",&
513 fouriername_t(:ilen(fouriername_t))
514 write (iout,*) "Torsional parameter file : ",&
515 torname_t(:ilen(torname_t))
516 write (iout,*) "Double torsional parameter file : ",&
517 tordname_t(:ilen(tordname_t))
518 write (iout,*) "Backbone-rotamer parameter file : ",&
519 sccorname(:ilen(sccorname))
520 write (iout,*) "Bond & inertia constant file : ",&
521 bondname_t(:ilen(bondname_t))
522 write (iout,*) "Bending parameter file : ",&
523 thetname_t(:ilen(thetname_t))
524 write (iout,*) "Rotamer parameter file : ",&
525 rotname_t(:ilen(rotname_t))
528 ! Read the virtual-bond parameters, masses, and moments of inertia
529 ! and Stokes' radii of the peptide group and side chains
531 allocate(dsc(ntyp1)) !(ntyp1)
532 allocate(dsc_inv(ntyp1)) !(ntyp1)
533 allocate(nbondterm(ntyp)) !(ntyp)
534 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
535 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
536 !el allocate(msc(ntyp+1)) !(ntyp+1)
537 !el allocate(isc(ntyp+1)) !(ntyp+1)
538 !el allocate(restok(ntyp+1)) !(ntyp+1)
539 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
542 read (ibond,*) vbldp0,akp
545 read (ibond,*) vbldsc0(1,i),aksc(1,i)
546 dsc(i) = vbldsc0(1,i)
550 dsc_inv(i)=1.0D0/dsc(i)
554 read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
556 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
558 dsc(i) = vbldsc0(1,i)
562 dsc_inv(i)=1.0D0/dsc(i)
567 write(iout,'(/a/)')"Force constants virtual bonds:"
568 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
570 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
572 write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
573 vbldsc0(1,i),aksc(1,i),abond0(1,i)
575 write (iout,'(13x,3f10.5)') &
576 vbldsc0(j,i),aksc(j,i),abond0(j,i)
580 if (.not. allocated(msc)) allocate(msc(ntyp1,5))
581 if (.not. allocated(restok)) allocate(restok(ntyp1,5))
584 read(iion,*) msc(i,5),restok(i,5)
585 print *,msc(i,5),restok(i,5)
589 read (iion,*) ncatprotparm
590 allocate(catprm(ncatprotparm,4))
592 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
597 !----------------------------------------------------
598 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
599 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
600 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
601 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
602 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
603 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
609 athet(j,i,ichir1,ichir2)=0.0D0
610 bthet(j,i,ichir1,ichir2)=0.0D0
624 !elwrite(iout,*) "parmread kontrol"
628 ! Read the parameters of the probability distribution/energy expression
629 ! of the virtual-bond valence angles theta
632 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
633 (bthet(j,i,1,1),j=1,2)
634 read (ithep,*) (polthet(j,i),j=0,3)
635 !elwrite(iout,*) "parmread kontrol in cryst_theta"
636 read (ithep,*) (gthet(j,i),j=1,3)
637 !elwrite(iout,*) "parmread kontrol in cryst_theta"
638 read (ithep,*) theta0(i),sig0(i),sigc0(i)
640 !elwrite(iout,*) "parmread kontrol in cryst_theta"
642 !elwrite(iout,*) "parmread kontrol in cryst_theta"
644 athet(1,i,1,-1)=athet(1,i,1,1)
645 athet(2,i,1,-1)=athet(2,i,1,1)
646 bthet(1,i,1,-1)=-bthet(1,i,1,1)
647 bthet(2,i,1,-1)=-bthet(2,i,1,1)
648 athet(1,i,-1,1)=-athet(1,i,1,1)
649 athet(2,i,-1,1)=-athet(2,i,1,1)
650 bthet(1,i,-1,1)=bthet(1,i,1,1)
651 bthet(2,i,-1,1)=bthet(2,i,1,1)
653 !elwrite(iout,*) "parmread kontrol in cryst_theta"
656 athet(1,i,-1,-1)=athet(1,-i,1,1)
657 athet(2,i,-1,-1)=-athet(2,-i,1,1)
658 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
659 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
660 athet(1,i,-1,1)=athet(1,-i,1,1)
661 athet(2,i,-1,1)=-athet(2,-i,1,1)
662 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
663 bthet(2,i,-1,1)=bthet(2,-i,1,1)
664 athet(1,i,1,-1)=-athet(1,-i,1,1)
665 athet(2,i,1,-1)=athet(2,-i,1,1)
666 bthet(1,i,1,-1)=bthet(1,-i,1,1)
667 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
672 polthet(j,i)=polthet(j,-i)
675 gthet(j,i)=gthet(j,-i)
678 !elwrite(iout,*) "parmread kontrol in cryst_theta"
680 !elwrite(iout,*) "parmread kontrol in cryst_theta"
683 ! & 'Parameters of the virtual-bond valence angles:'
684 ! write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
685 ! & ' ATHETA0 ',' A1 ',' A2 ',
688 ! write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
689 ! & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
691 ! write (iout,'(/a/9x,5a/79(1h-))')
692 ! & 'Parameters of the expression for sigma(theta_c):',
693 ! & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
694 ! & ' ALPH3 ',' SIGMA0C '
696 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
697 ! & (polthet(j,i),j=0,3),sigc0(i)
699 ! write (iout,'(/a/9x,5a/79(1h-))')
700 ! & 'Parameters of the second gaussian:',
701 ! & ' THETA0 ',' SIGMA0 ',' G1 ',
704 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
705 ! & sig0(i),(gthet(j,i),j=1,3)
708 'Parameters of the virtual-bond valence angles:'
709 write (iout,'(/a/9x,5a/79(1h-))') &
710 'Coefficients of expansion',&
711 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
712 ' b1*10^1 ',' b2*10^1 '
714 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
715 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
716 (10*bthet(j,i,1,1),j=1,2)
718 write (iout,'(/a/9x,5a/79(1h-))') &
719 'Parameters of the expression for sigma(theta_c):',&
720 ' alpha0 ',' alph1 ',' alph2 ',&
721 ' alhp3 ',' sigma0c '
723 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
724 (polthet(j,i),j=0,3),sigc0(i)
726 write (iout,'(/a/9x,5a/79(1h-))') &
727 'Parameters of the second gaussian:',&
728 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
731 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
732 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
737 ! Read the parameters of Utheta determined from ab initio surfaces
738 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
740 ! write (iout,*) "tu dochodze"
741 read (ithep,*) nthetyp,ntheterm,ntheterm2,&
742 ntheterm3,nsingle,ndouble
743 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
745 !----------------------------------------------------
746 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
747 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
748 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
749 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
750 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
751 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
752 !(maxtheterm,-maxthetyp1:maxthetyp1,&
753 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
754 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
755 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
756 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
757 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
758 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
759 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
760 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
761 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
762 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
763 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
764 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
765 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
766 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
767 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
768 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
769 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
772 read (ithep,*) (ithetyp(i),i=1,ntyp1)
774 ithetyp(i)=-ithetyp(-i)
776 ! write (iout,*) "tu dochodze"
777 aa0thet(:,:,:,:)=0.0d0
778 aathet(:,:,:,:,:)=0.0d0
779 bbthet(:,:,:,:,:,:)=0.0d0
780 ccthet(:,:,:,:,:,:)=0.0d0
781 ddthet(:,:,:,:,:,:)=0.0d0
782 eethet(:,:,:,:,:,:)=0.0d0
783 ffthet(:,:,:,:,:,:,:)=0.0d0
784 ggthet(:,:,:,:,:,:,:)=0.0d0
788 do j=-nthetyp,nthetyp
789 do k=-nthetyp,nthetyp
790 read (ithep,'(6a)') res1
791 read (ithep,*) aa0thet(i,j,k,iblock)
792 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
794 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
795 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
796 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
797 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
800 (((ffthet(llll,lll,ll,i,j,k,iblock),&
801 ffthet(lll,llll,ll,i,j,k,iblock),&
802 ggthet(llll,lll,ll,i,j,k,iblock),&
803 ggthet(lll,llll,ll,i,j,k,iblock),&
804 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
809 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
810 ! coefficients of theta-and-gamma-dependent terms are zero.
815 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
816 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
818 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
819 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
822 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
824 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
827 ! Substitution for D aminoacids from symmetry.
830 do j=-nthetyp,nthetyp
831 do k=-nthetyp,nthetyp
832 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
834 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
838 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
839 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
840 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
841 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
847 ffthet(llll,lll,ll,i,j,k,iblock)= &
848 ffthet(llll,lll,ll,-i,-j,-k,iblock)
849 ffthet(lll,llll,ll,i,j,k,iblock)= &
850 ffthet(lll,llll,ll,-i,-j,-k,iblock)
851 ggthet(llll,lll,ll,i,j,k,iblock)= &
852 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
853 ggthet(lll,llll,ll,i,j,k,iblock)= &
854 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
864 ! Control printout of the coefficients of virtual-bond-angle potentials
868 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
872 write (iout,'(//4a)') &
873 'Type ',onelett(i),onelett(j),onelett(k)
874 write (iout,'(//a,10x,a)') " l","a[l]"
875 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
876 write (iout,'(i2,1pe15.5)') &
877 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
879 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
880 "b",l,"c",l,"d",l,"e",l
882 write (iout,'(i2,4(1pe15.5))') m,&
883 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
884 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
888 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
889 "f+",l,"f-",l,"g+",l,"g-",l
892 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
893 ffthet(n,m,l,i,j,k,iblock),&
894 ffthet(m,n,l,i,j,k,iblock),&
895 ggthet(n,m,l,i,j,k,iblock),&
896 ggthet(m,n,l,i,j,k,iblock)
907 !-------------------------------------------
908 allocate(nlob(ntyp1)) !(ntyp1)
909 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
910 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
911 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
929 gaussc(l,k,j,i)=0.0D0
937 ! Read the parameters of the probability distribution/energy expression
938 ! of the side chains.
941 !c write (iout,*) "tu dochodze",i
942 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
946 dsc_inv(i)=1.0D0/dsc(i)
957 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
958 censc(1,1,-i)=censc(1,1,i)
959 censc(2,1,-i)=censc(2,1,i)
960 censc(3,1,-i)=-censc(3,1,i)
962 read (irotam,*) bsc(j,i)
963 read (irotam,*) (censc(k,j,i),k=1,3),&
964 ((blower(k,l,j),l=1,k),k=1,3)
965 censc(1,j,-i)=censc(1,j,i)
966 censc(2,j,-i)=censc(2,j,i)
967 censc(3,j,-i)=-censc(3,j,i)
968 ! BSC is amplitude of Gaussian
975 akl=akl+blower(k,m,j)*blower(l,m,j)
979 if (((k.eq.3).and.(l.ne.3)) &
980 .or.((l.eq.3).and.(k.ne.3))) then
981 gaussc(k,l,j,-i)=-akl
982 gaussc(l,k,j,-i)=-akl
994 write (iout,'(/a)') 'Parameters of side-chain local geometry'
998 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
999 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1000 ! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1001 ! write (iout,'(a,f10.4,4(16x,f10.4))')
1002 ! & 'Center ',(bsc(j,i),j=1,nlobi)
1003 ! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
1004 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1005 'log h',(bsc(j,i),j=1,nlobi)
1006 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1007 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1008 ! write (iout,'(a)')
1014 ! blower(k,l,j)=gaussc(ind,j,i)
1019 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1020 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1027 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1028 ! added by Urszula Kozlowska 07/11/2007
1030 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1038 read(irotam,*) sc_parmin(j,i)
1046 ! Read torsional parameters in old format
1048 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1050 read (itorp,*) ntortyp,nterm_old
1051 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1052 read (itorp,*) (itortyp(i),i=1,ntyp)
1054 !el from energy module--------
1055 allocate(v1(nterm_old,ntortyp,ntortyp))
1056 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1057 !el---------------------------
1063 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1069 write (iout,'(/a/)') 'Torsional constants:'
1072 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1073 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1081 ! Read torsional parameters
1083 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1085 read (itorp,*) ntortyp
1086 read (itorp,*) (itortyp(i),i=1,ntyp)
1087 write (iout,*) 'ntortyp',ntortyp
1089 !el from energy module---------
1090 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1091 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1093 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1094 allocate(vlor2(maxlor,ntortyp,ntortyp))
1095 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1096 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1098 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1099 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1100 !el---------------------------
1102 do i=-ntortyp,ntortyp
1103 do j=-ntortyp,ntortyp
1109 !el---------------------------
1113 itortyp(i)=-itortyp(-i)
1115 ! write (iout,*) 'ntortyp',ntortyp
1117 do j=-ntortyp+1,ntortyp-1
1118 read (itorp,*) nterm(i,j,iblock),&
1120 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1121 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1124 do k=1,nterm(i,j,iblock)
1125 read (itorp,*) kk,v1(k,i,j,iblock),&
1127 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1128 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1129 v0ij=v0ij+si*v1(k,i,j,iblock)
1132 do k=1,nlor(i,j,iblock)
1133 read (itorp,*) kk,vlor1(k,i,j),&
1134 vlor2(k,i,j),vlor3(k,i,j)
1135 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1138 v0(-i,-j,iblock)=v0ij
1145 write (iout,'(/a/)') 'Torsional constants:'
1148 write (iout,*) 'ityp',i,' jtyp',j
1149 write (iout,*) 'Fourier constants'
1150 do k=1,nterm(i,j,iblock)
1151 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1154 write (iout,*) 'Lorenz constants'
1155 do k=1,nlor(i,j,iblock)
1156 write (iout,'(3(1pe15.5))') &
1157 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1164 ! 6/23/01 Read parameters for double torsionals
1166 !el from energy module------------
1167 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1168 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1169 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1170 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1171 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1172 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1173 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1174 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1175 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1176 !---------------------------------
1180 do j=-ntortyp+1,ntortyp-1
1181 do k=-ntortyp+1,ntortyp-1
1182 read (itordp,'(3a1)') t1,t2,t3
1183 ! write (iout,*) "OK onelett",
1186 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1187 .or. t3.ne.toronelet(k)) then
1188 write (iout,*) "Error in double torsional parameter file",&
1191 call MPI_Finalize(Ierror)
1193 stop "Error in double torsional parameter file"
1195 read (itordp,*) ntermd_1(i,j,k,iblock),&
1196 ntermd_2(i,j,k,iblock)
1197 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1198 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1199 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1200 ntermd_1(i,j,k,iblock))
1201 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1202 ntermd_1(i,j,k,iblock))
1203 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1204 ntermd_1(i,j,k,iblock))
1205 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1206 ntermd_1(i,j,k,iblock))
1207 ! Martix of D parameters for one dimesional foureir series
1208 do l=1,ntermd_1(i,j,k,iblock)
1209 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1210 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1211 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1212 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1213 ! write(iout,*) "whcodze" ,
1214 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1216 read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1217 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1218 v2s(m,l,i,j,k,iblock),&
1219 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1220 ! Martix of D parameters for two dimesional fourier series
1221 do l=1,ntermd_2(i,j,k,iblock)
1223 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1224 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1225 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1226 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1235 write (iout,*) 'Constants for double torsionals'
1238 do j=-ntortyp+1,ntortyp-1
1239 do k=-ntortyp+1,ntortyp-1
1240 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1241 ' nsingle',ntermd_1(i,j,k,iblock),&
1242 ' ndouble',ntermd_2(i,j,k,iblock)
1244 write (iout,*) 'Single angles:'
1245 do l=1,ntermd_1(i,j,k,iblock)
1246 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1247 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1248 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1249 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1252 write (iout,*) 'Pairs of angles:'
1253 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1254 do l=1,ntermd_2(i,j,k,iblock)
1255 write (iout,'(i5,20f10.5)') &
1256 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1259 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1260 do l=1,ntermd_2(i,j,k,iblock)
1261 write (iout,'(i5,20f10.5)') &
1262 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1263 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1272 !elwrite(iout,*) "parmread kontrol sc-bb"
1273 ! Read of Side-chain backbone correlation parameters
1274 ! Modified 11 May 2012 by Adasko
1277 read (isccor,*) nsccortyp
1280 !c maxinter is maximum interaction sites
1281 !write(iout,*)"maxterm_sccor",maxterm_sccor
1282 !el from module energy-------------
1283 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1284 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1285 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1286 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1287 !-----------------------------------
1288 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1289 !-----------------------------------
1290 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1291 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1292 -nsccortyp:nsccortyp))
1293 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1294 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1295 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1296 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1297 !-----------------------------------
1298 do i=-nsccortyp,nsccortyp
1299 do j=-nsccortyp,nsccortyp
1303 !-----------------------------------
1305 read (isccor,*) (isccortyp(i),i=1,ntyp)
1307 isccortyp(i)=-isccortyp(-i)
1309 iscprol=isccortyp(20)
1310 ! write (iout,*) 'ntortyp',ntortyp
1312 !c maxinter is maximum interaction sites
1317 nterm_sccor(i,j),nlor_sccor(i,j)
1323 nterm_sccor(-i,j)=nterm_sccor(i,j)
1324 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1325 nterm_sccor(i,-j)=nterm_sccor(i,j)
1326 do k=1,nterm_sccor(i,j)
1327 read (isccor,*) kk,v1sccor(k,l,i,j),&
1329 if (j.eq.iscprol) then
1330 if (i.eq.isccortyp(10)) then
1331 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1332 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1334 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1335 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1336 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1337 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1338 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1339 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1340 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1341 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1344 if (i.eq.isccortyp(10)) then
1345 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1346 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1348 if (j.eq.isccortyp(10)) then
1349 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1350 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1352 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1353 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1354 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1355 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1356 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1357 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1361 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1362 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1363 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1364 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1367 do k=1,nlor_sccor(i,j)
1368 read (isccor,*) kk,vlor1sccor(k,i,j),&
1369 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1370 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1371 (1+vlor3sccor(k,i,j)**2)
1373 v0sccor(l,i,j)=v0ijsccor
1374 v0sccor(l,-i,j)=v0ijsccor1
1375 v0sccor(l,i,-j)=v0ijsccor2
1376 v0sccor(l,-i,-j)=v0ijsccor3
1382 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1385 write (iout,*) 'ityp',i,' jtyp',j
1386 write (iout,*) 'Fourier constants'
1387 do k=1,nterm_sccor(i,j)
1388 write (iout,'(2(1pe15.5))') &
1389 (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
1391 write (iout,*) 'Lorenz constants'
1392 do k=1,nlor_sccor(i,j)
1393 write (iout,'(3(1pe15.5))') &
1394 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1400 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1401 ! interaction energy of the Gly, Ala, and Pro prototypes.
1403 read (ifourier,*) nloctyp
1404 !el write(iout,*)"nloctyp",nloctyp
1405 !el from module energy-------
1406 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1407 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1408 allocate(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1409 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1410 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1411 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1412 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1413 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1415 do ii=-nloctyp-1,nloctyp+1
1423 ctilde(j,i,ii)=0.0d0
1424 dtilde(j,i,ii)=0.0d0
1428 !--------------------------------
1429 allocate(b(13,0:nloctyp))
1433 read (ifourier,*) (b(ii,i),ii=1,13)
1435 write (iout,*) 'Type',i
1436 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1444 B1tilde(1,i) = b(3,i)
1445 B1tilde(2,i) =-b(5,i)
1446 B1tilde(1,-i) =-b(3,i)
1447 B1tilde(2,-i) =b(5,i)
1448 ! b1tilde(1,i)=0.0d0
1449 ! b1tilde(2,i)=0.0d0
1469 Ctilde(1,1,i)=b(7,i)
1470 Ctilde(1,2,i)=b(9,i)
1471 Ctilde(2,1,i)=-b(9,i)
1472 Ctilde(2,2,i)=b(7,i)
1473 Ctilde(1,1,-i)=b(7,i)
1474 Ctilde(1,2,-i)=-b(9,i)
1475 Ctilde(2,1,-i)=b(9,i)
1476 Ctilde(2,2,-i)=b(7,i)
1478 ! Ctilde(1,1,i)=0.0d0
1479 ! Ctilde(1,2,i)=0.0d0
1480 ! Ctilde(2,1,i)=0.0d0
1481 ! Ctilde(2,2,i)=0.0d0
1494 Dtilde(1,1,i)=b(6,i)
1495 Dtilde(1,2,i)=b(8,i)
1496 Dtilde(2,1,i)=-b(8,i)
1497 Dtilde(2,2,i)=b(6,i)
1498 Dtilde(1,1,-i)=b(6,i)
1499 Dtilde(1,2,-i)=-b(8,i)
1500 Dtilde(2,1,-i)=b(8,i)
1501 Dtilde(2,2,-i)=b(6,i)
1503 ! Dtilde(1,1,i)=0.0d0
1504 ! Dtilde(1,2,i)=0.0d0
1505 ! Dtilde(2,1,i)=0.0d0
1506 ! Dtilde(2,2,i)=0.0d0
1507 EE(1,1,i)= b(10,i)+b(11,i)
1508 EE(2,2,i)=-b(10,i)+b(11,i)
1509 EE(2,1,i)= b(12,i)-b(13,i)
1510 EE(1,2,i)= b(12,i)+b(13,i)
1511 EE(1,1,-i)= b(10,i)+b(11,i)
1512 EE(2,2,-i)=-b(10,i)+b(11,i)
1513 EE(2,1,-i)=-b(12,i)+b(13,i)
1514 EE(1,2,-i)=-b(12,i)-b(13,i)
1520 ! ee(2,1,i)=ee(1,2,i)
1525 write (iout,*) 'Type',i
1527 ! write (iout,'(f10.5)') B1(:,i)
1528 write(iout,*) B1(1,i),B1(2,i)
1530 ! write (iout,'(f10.5)') B2(:,i)
1531 write(iout,*) B2(1,i),B2(2,i)
1534 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1538 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1542 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1547 ! Read electrostatic-interaction parameters
1550 write (iout,'(/a)') 'Electrostatic interaction constants:'
1551 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1552 'IT','JT','APP','BPP','AEL6','AEL3'
1554 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1555 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1556 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1557 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1562 app (i,j)=epp(i,j)*rri*rri
1563 bpp (i,j)=-2.0D0*epp(i,j)*rri
1564 ael6(i,j)=elpp6(i,j)*4.2D0**6
1565 ael3(i,j)=elpp3(i,j)*4.2D0**3
1566 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1571 ! Read side-chain interaction parameters.
1573 !el from module energy - COMMON.INTERACT-------
1574 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1575 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1576 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1577 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1578 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1579 allocate(epslip(ntyp,ntyp))
1590 !--------------------------------
1592 read (isidep,*) ipot,expon
1593 !el if (ipot.lt.1 .or. ipot.gt.5) then
1594 ! write (iout,'(2a)') 'Error while reading SC interaction',&
1595 ! 'potential file - unknown potential type.'
1599 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1600 ', exponents are ',expon,2*expon
1601 ! goto (10,20,30,30,40) ipot
1603 !----------------------- LJ potential ---------------------------------
1605 ! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1606 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1608 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1609 write (iout,'(a/)') 'The epsilon array:'
1610 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1611 write (iout,'(/a)') 'One-body parameters:'
1612 write (iout,'(a,4x,a)') 'residue','sigma'
1613 write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
1616 !----------------------- LJK potential --------------------------------
1618 ! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1619 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1620 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1622 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1623 write (iout,'(a/)') 'The epsilon array:'
1624 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1625 write (iout,'(/a)') 'One-body parameters:'
1626 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1627 write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
1631 !---------------------- GB or BP potential -----------------------------
1635 read (isidep,*)(eps(i,j),j=i,ntyp)
1637 read (isidep,*)(sigma0(i),i=1,ntyp)
1638 read (isidep,*)(sigii(i),i=1,ntyp)
1639 read (isidep,*)(chip(i),i=1,ntyp)
1640 read (isidep,*)(alp(i),i=1,ntyp)
1642 read (isidep,*)(epslip(i,j),j=i,ntyp)
1644 ! For the GB potential convert sigma'**2 into chi'
1647 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1651 write (iout,'(/a/)') 'Parameters of the BP potential:'
1652 write (iout,'(a/)') 'The epsilon array:'
1653 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1654 write (iout,'(/a)') 'One-body parameters:'
1655 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
1657 write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
1658 chip(i),alp(i),i=1,ntyp)
1661 !--------------------- GBV potential -----------------------------------
1663 ! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1664 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1665 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
1666 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1668 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1669 write (iout,'(a/)') 'The epsilon array:'
1670 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1671 write (iout,'(/a)') 'One-body parameters:'
1672 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
1673 's||/s_|_^2',' chip ',' alph '
1674 write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
1675 sigii(i),chip(i),alp(i),i=1,ntyp)
1678 write (iout,'(2a)') 'Error while reading SC interaction',&
1679 'potential file - unknown potential type.'
1685 !-----------------------------------------------------------------------
1686 ! Calculate the "working" parameters of SC interactions.
1688 !el from module energy - COMMON.INTERACT-------
1689 ! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
1690 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
1691 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
1692 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
1704 !--------------------------------
1709 epslip(i,j)=epslip(j,i)
1714 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1715 sigma(j,i)=sigma(i,j)
1716 rs0(i,j)=dwa16*sigma(i,j)
1720 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
1721 'Working parameters of the SC interactions:',&
1722 ' a ',' b ',' augm ',' sigma ',' r0 ',&
1727 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1736 sigeps=dsign(1.0D0,epsij)
1738 aa_aq(i,j)=epsij*rrij*rrij
1739 bb_aq(i,j)=-sigeps*epsij*rrij
1740 aa_aq(j,i)=aa_aq(i,j)
1741 bb_aq(j,i)=bb_aq(i,j)
1742 epsijlip=epslip(i,j)
1743 sigeps=dsign(1.0D0,epsijlip)
1744 epsijlip=dabs(epsijlip)
1745 aa_lip(i,j)=epsijlip*rrij*rrij
1746 bb_lip(i,j)=-sigeps*epsijlip*rrij
1747 aa_lip(j,i)=aa_lip(i,j)
1748 bb_lip(j,i)=bb_lip(i,j)
1750 sigt1sq=sigma0(i)**2
1751 sigt2sq=sigma0(j)**2
1754 ratsig1=sigt2sq/sigt1sq
1755 ratsig2=1.0D0/ratsig1
1756 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1757 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1758 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1762 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1763 sigmaii(i,j)=rsum_max
1764 sigmaii(j,i)=rsum_max
1766 ! sigmaii(i,j)=r0(i,j)
1767 ! sigmaii(j,i)=r0(i,j)
1769 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1770 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1771 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1772 augm(i,j)=epsij*r_augm**(2*expon)
1773 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1780 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
1781 restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
1782 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1787 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
1795 ! Define the SC-p interaction constants
1805 !elwrite(iout,*) "parmread kontrol before oldscp"
1807 ! Define the SC-p interaction constants
1811 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1813 ! aad(i,1)=0.3D0*4.0D0**12
1814 ! Following line for constants currently implemented
1815 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
1816 aad(i,1)=1.5D0*4.0D0**12
1817 ! aad(i,1)=0.17D0*5.6D0**12
1819 ! "Soft" SC-p repulsion
1821 ! Following line for constants currently implemented
1822 ! aad(i,1)=0.3D0*4.0D0**6
1823 ! "Hard" SC-p repulsion
1824 bad(i,1)=3.0D0*4.0D0**6
1825 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1834 ! 8/9/01 Read the SC-p interaction constants from file
1837 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1840 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1841 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1842 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1843 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1847 write (iout,*) "Parameters of SC-p interactions:"
1849 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
1850 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1855 ! Define the constants of the disulfide bridge
1859 ! Old arbitrary potential - commented out.
1864 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1865 ! energy surface of diethyl disulfide.
1866 ! A. Liwo and U. Kozlowska, 11/24/03
1877 write (iout,'(/a)') "Disulfide bridge parameters:"
1878 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1879 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1880 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1881 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
1885 end subroutine parmread
1887 !-----------------------------------------------------------------------------
1889 !-----------------------------------------------------------------------------
1890 subroutine mygetenv(string,var)
1894 ! This subroutine passes the environmental variables to FORTRAN program.
1895 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
1896 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
1897 ! reads the environmental variables from $HOME/.env
1899 ! Usage: As for the standard FORTRAN GETENV subroutine.
1901 ! Purpose: some versions/installations of MPI do not transfer the environmental
1902 ! variables to slave processors, if these variables are set in the shell script
1903 ! from which mpirun is called.
1912 character*(*) :: string,var
1913 #if defined(MYGETENV) && defined(MPI)
1914 ! include "DIMENSIONS.ZSCOPT"
1916 ! include "COMMON.MPI"
1917 !el character*360 ucase
1919 character(len=360) :: string1(360),karta
1920 character(len=240) :: home
1923 call getenv("HOME",home)
1924 open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
1926 read (99,end=111,err=111,'(a)') karta
1930 call split_string(karta,string1,80,n)
1931 if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
1932 string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
1934 print *,"Processor",me,": ",var(:ilen(var)),&
1935 " assigned to ",string(:ilen(string))
1940 111 print *,"Environment variable ",string(:ilen(string))," not set."
1943 112 print *,"Error opening environment file!"
1945 call getenv(string,var)
1948 end subroutine mygetenv
1949 !-----------------------------------------------------------------------------
1951 !-----------------------------------------------------------------------------
1952 subroutine read_general_data(*)
1954 use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions
1955 use energy_data, only:distchainmax
1956 use geometry_data, only:boxxsize,boxysize,boxzsize
1958 ! include "DIMENSIONS"
1959 ! include "DIMENSIONS.ZSCOPT"
1960 ! include "DIMENSIONS.FREE"
1961 ! include "COMMON.TORSION"
1962 ! include "COMMON.INTERACT"
1963 ! include "COMMON.IOUNITS"
1964 ! include "COMMON.TIME1"
1965 ! include "COMMON.PROT"
1966 ! include "COMMON.PROTFILES"
1967 ! include "COMMON.CHAIN"
1968 ! include "COMMON.NAMES"
1969 ! include "COMMON.FFIELD"
1970 ! include "COMMON.ENEPS"
1971 ! include "COMMON.WEIGHTS"
1972 ! include "COMMON.FREE"
1973 ! include "COMMON.CONTROL"
1974 ! include "COMMON.ENERGIES"
1975 character(len=800) :: controlcard
1976 integer :: i,j,k,ii,n_ene_found
1977 integer :: ind,itype1,itype2,itypf,itypsc,itypp
1980 !el character*16 ucase
1981 character(len=16) :: key
1983 call card_concat(controlcard,.true.)
1984 call readi(controlcard,"N_ENE",n_eneW,max_eneW)
1985 if (n_eneW.gt.max_eneW) then
1986 write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
1990 call readi(controlcard,"NPARMSET",nparmset,1)
1991 !elwrite(iout,*)"in read_gen data"
1992 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
1993 call readi(controlcard,"IPARMPRINT",iparmprint,1)
1994 write (iout,*) "PARMPRINT",iparmprint
1995 if (nparmset.gt.max_parm) then
1996 write (iout,*) "Error: parameter out of range: NPARMSET",&
2000 !elwrite(iout,*)"in read_gen data"
2001 call readi(controlcard,"MAXIT",maxit,5000)
2002 call reada(controlcard,"FIMIN",fimin,1.0d-3)
2003 call readi(controlcard,"ENSEMBLES",ensembles,0)
2004 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
2005 write (iout,*) "Number of energy parameter sets",nparmset
2006 allocate(isampl(nparmset))
2007 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
2008 write (iout,*) "MaxSlice",MaxSlice
2009 call readi(controlcard,"NSLICE",nslice,1)
2010 !elwrite(iout,*)"in read_gen data"
2012 if (nslice.gt.MaxSlice) then
2013 write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
2017 write (iout,*) "Frequency of storing conformations",&
2018 (isampl(i),i=1,nparmset)
2019 write (iout,*) "Maxit",maxit," Fimin",fimin
2020 call readi(controlcard,"NQ",nQ,1)
2021 if (nQ.gt.MaxQ) then
2022 write (iout,*) "Error: parameter out of range: NQ",nq,&
2027 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
2028 call reada(controlcard,"DELTA",delta,1.0d-2)
2029 call readi(controlcard,"EINICHECK",einicheck,2)
2030 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
2031 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
2032 call readi(controlcard,"RESCALE",rescale_modeW,1)
2033 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2034 call reada(controlcard,'BOXY',boxysize,100.0d0)
2035 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2036 ions=index(controlcard,"IONS").gt.0
2037 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2038 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2039 write(iout,*) "R_CUT_ELE=",r_cut_ele
2040 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
2041 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
2042 call readi(controlcard,'SYM',symetr,1)
2043 write (iout,*) "DISTCHAINMAX",distchainmax
2044 write (iout,*) "delta",delta
2045 write (iout,*) "einicheck",einicheck
2046 write (iout,*) "rescale_mode",rescale_modeW
2048 bxfile=index(controlcard,"BXFILE").gt.0
2049 cxfile=index(controlcard,"CXFILE").gt.0
2050 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
2052 histfile=index(controlcard,"HISTFILE").gt.0
2053 histout=index(controlcard,"HISTOUT").gt.0
2054 entfile=index(controlcard,"ENTFILE").gt.0
2055 zscfile=index(controlcard,"ZSCFILE").gt.0
2056 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
2057 write (iout,*) "with_dihed_constr ",with_dihed_constr
2058 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2060 end subroutine read_general_data
2061 !------------------------------------------------------------------------------
2062 subroutine read_efree(*)
2064 ! Read molecular data
2067 ! include 'DIMENSIONS'
2068 ! include 'DIMENSIONS.ZSCOPT'
2069 ! include 'DIMENSIONS.COMPAR'
2070 ! include 'DIMENSIONS.FREE'
2071 ! include 'COMMON.IOUNITS'
2072 ! include 'COMMON.TIME1'
2073 ! include 'COMMON.SBRIDGE'
2074 ! include 'COMMON.CONTROL'
2075 ! include 'COMMON.CHAIN'
2076 ! include 'COMMON.HEADER'
2077 ! include 'COMMON.GEO'
2078 ! include 'COMMON.FREE'
2079 character(len=320) :: controlcard !,ucase
2080 integer :: iparm,ib,i,j,npars
2090 ! call alloc_wham_arrays
2091 ! allocate(nT_h(nParmSet))
2092 ! allocate(replica(nParmSet))
2093 ! allocate(umbrella(nParmSet))
2094 ! allocate(read_iset(nParmSet))
2095 ! allocate(nT_h(nParmSet))
2099 call card_concat(controlcard,.true.)
2100 call readi(controlcard,'NT',nT_h(iparm),1)
2101 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
2103 if (nT_h(iparm).gt.MaxT_h) then
2104 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
2108 replica(iparm)=index(controlcard,"REPLICA").gt.0
2109 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
2110 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
2111 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
2112 replica(iparm)," umbrella ",umbrella(iparm),&
2113 " read_iset",read_iset(iparm)
2116 call card_concat(controlcard,.true.)
2117 call readi(controlcard,'NR',nR(ib,iparm),1)
2118 if (umbrella(iparm)) then
2121 nRR(ib,iparm)=nR(ib,iparm)
2123 if (nR(ib,iparm).gt.MaxR) then
2124 write (iout,*) "Error: parameter out of range: NR",&
2128 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
2129 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
2130 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
2133 call card_concat(controlcard,.true.)
2134 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
2136 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
2141 write (iout,*) "ib",ib," beta_h",&
2142 1.0d0/(0.001987*beta_h(ib,iparm))
2143 write (iout,*) "nR",nR(ib,iparm)
2144 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
2146 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
2147 "q0",(q0(j,i,ib,iparm),j=1,nQ)
2159 nR(ib,iparm)=nR(ib,1)
2160 if (umbrella(iparm)) then
2163 nRR(ib,iparm)=nR(ib,1)
2165 beta_h(ib,iparm)=beta_h(ib,1)
2167 f(i,ib,iparm)=f(i,ib,1)
2169 KH(j,i,ib,iparm)=KH(j,i,ib,1)
2170 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
2173 replica(iparm)=replica(1)
2174 umbrella(iparm)=umbrella(1)
2175 read_iset(iparm)=read_iset(1)
2182 end subroutine read_efree
2183 !-----------------------------------------------------------------------------
2184 subroutine read_protein_data(*)
2186 ! include "DIMENSIONS"
2187 ! include "DIMENSIONS.ZSCOPT"
2188 ! include "DIMENSIONS.FREE"
2192 integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
2193 ! include "COMMON.MPI"
2195 ! include "COMMON.CHAIN"
2196 ! include "COMMON.IOUNITS"
2197 ! include "COMMON.PROT"
2198 ! include "COMMON.PROTFILES"
2199 ! include "COMMON.NAMES"
2200 ! include "COMMON.FREE"
2201 ! include "COMMON.OBCINKA"
2202 character(len=64) :: nazwa
2203 character(len=16000) :: controlcard
2204 integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
2205 !el external ilen,iroof
2214 ! Read names of files with conformation data.
2215 if (replica(iparm)) then
2221 do ii=1,nRR(ib,iparm)
2222 write (iout,*) "Parameter set",iparm," temperature",ib,&
2225 call card_concat(controlcard,.true.)
2226 write (iout,*) controlcard(:ilen(controlcard))
2227 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
2228 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
2229 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
2230 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
2231 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
2232 maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
2233 call reada(controlcard,"TIME_START",&
2234 time_start_collect(ii,ib,iparm),0.0d0)
2235 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
2237 write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
2238 " rec_end",rec_end(ii,ib,iparm)
2239 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
2240 " time_end",time_end_collect(ii,ib,iparm)
2242 if (replica(iparm)) then
2243 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
2244 write (iout,*) "Number of trajectories",totraj(ii,iparm)
2247 if (nfile_bin(ii,ib,iparm).lt.2 &
2248 .and. nfile_asc(ii,ib,iparm).eq.0 &
2249 .and. nfile_cx(ii,ib,iparm).eq.0) then
2250 write (iout,*) "Error - no action specified!"
2253 if (nfile_bin(ii,ib,iparm).gt.0) then
2254 call card_concat(controlcard,.false.)
2255 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
2256 maxfile_prot,nfile_bin(ii,ib,iparm))
2258 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
2259 write(iout,*) (protfiles(i,1,ii,ib,iparm),&
2260 i=1,nfile_bin(ii,ib,iparm))
2263 if (nfile_asc(ii,ib,iparm).gt.0) then
2264 call card_concat(controlcard,.false.)
2265 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2266 maxfile_prot,nfile_asc(ii,ib,iparm))
2268 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
2269 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2270 i=1,nfile_asc(ii,ib,iparm))
2272 else if (nfile_cx(ii,ib,iparm).gt.0) then
2273 call card_concat(controlcard,.false.)
2274 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2275 maxfile_prot,nfile_cx(ii,ib,iparm))
2277 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
2278 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2279 i=1,nfile_cx(ii,ib,iparm))
2288 end subroutine read_protein_data
2289 !-------------------------------------------------------------------------------
2290 subroutine readsss(rekord,lancuch,wartosc,default)
2292 character*(*) :: rekord,lancuch,wartosc,default
2293 character(len=80) :: aux
2294 integer :: lenlan,lenrec,iread,ireade
2298 lenlan=ilen(lancuch)
2300 iread=index(rekord,lancuch(:lenlan)//"=")
2301 ! print *,"rekord",rekord," lancuch",lancuch
2302 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2303 if (iread.eq.0) then
2307 iread=iread+lenlan+1
2308 ! print *,"iread",iread
2309 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2310 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2312 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2314 ! print *,"iread",iread
2315 if (iread.gt.lenrec) then
2320 ! print *,"ireade",ireade
2321 do while (ireade.lt.lenrec .and. &
2322 .not.iblnk(rekord(ireade:ireade)))
2325 wartosc=rekord(iread:ireade)
2327 end subroutine readsss
2328 !----------------------------------------------------------------------------
2329 subroutine multreads(rekord,lancuch,tablica,dim,default)
2332 character*(*) rekord,lancuch,tablica(dim),default
2333 character(len=80) :: aux
2334 integer :: lenlan,lenrec,iread,ireade
2341 lenlan=ilen(lancuch)
2343 iread=index(rekord,lancuch(:lenlan)//"=")
2344 ! print *,"rekord",rekord," lancuch",lancuch
2345 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2346 if (iread.eq.0) return
2347 iread=iread+lenlan+1
2349 ! print *,"iread",iread
2350 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2351 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2353 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2355 ! print *,"iread",iread
2356 if (iread.gt.lenrec) return
2358 ! print *,"ireade",ireade
2359 do while (ireade.lt.lenrec .and. &
2360 .not.iblnk(rekord(ireade:ireade)))
2363 tablica(i)=rekord(iread:ireade)
2366 end subroutine multreads
2367 !----------------------------------------------------------------------------
2368 subroutine split_string(rekord,tablica,dim,nsub)
2370 integer :: dim,nsub,i,ii,ll,kk
2371 character*(*) tablica(dim)
2372 character*(*) rekord
2382 ! Find the start of term name
2384 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2387 ! Parse the name into TABLICA(i) until blank found
2388 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2390 tablica(i)(kk:kk)=rekord(ii:ii)
2393 if (kk.gt.0) nsub=nsub+1
2394 if (ii.gt.ll) return
2397 end subroutine split_string
2398 !--------------------------------------------------------------------------------
2400 !--------------------------------------------------------------------------------
2401 subroutine read_compar
2403 ! Read molecular data
2405 use conform_compar, only:alloc_compar_arrays
2406 use control_data, only:pdbref
2407 use geometry_data, only:deg2rad,rad2deg
2409 ! include 'DIMENSIONS'
2410 ! include 'DIMENSIONS.ZSCOPT'
2411 ! include 'DIMENSIONS.COMPAR'
2412 ! include 'DIMENSIONS.FREE'
2413 ! include 'COMMON.IOUNITS'
2414 ! include 'COMMON.TIME1'
2415 ! include 'COMMON.SBRIDGE'
2416 ! include 'COMMON.CONTROL'
2417 ! include 'COMMON.COMPAR'
2418 ! include 'COMMON.CHAIN'
2419 ! include 'COMMON.HEADER'
2420 ! include 'COMMON.GEO'
2421 ! include 'COMMON.FREE'
2422 character(len=320) :: controlcard !,ucase
2423 character(len=64) :: wfile
2427 !elwrite(iout,*)"jestesmy w read_compar"
2428 call card_concat(controlcard,.true.)
2429 pdbref=(index(controlcard,'PDBREF').gt.0)
2430 call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
2431 call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
2432 call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
2433 call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
2434 verbose = index(controlcard,"VERBOSE").gt.0
2435 lgrp=index(controlcard,"STATIN").gt.0
2436 lgrp_out=index(controlcard,"STATOUT").gt.0
2437 merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
2438 binary = index(controlcard,"BINARY").gt.0
2439 rmscut_base_up=rmscut_base_up/50
2440 rmscut_base_low=rmscut_base_low/50
2441 call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
2442 call readi(controlcard,'NLEVEL',nlevel,1)
2443 if (nlevel.lt.0) then
2445 call alloc_compar_arrays(maxfrag,1)
2448 allocate(nfrag(nlevel))
2450 ! Read the data pertaining to elementary fragments (level 1)
2451 call readi(controlcard,'NFRAG',nfrag(1),0)
2452 write(iout,*)"nfrag(1)",nfrag(1)
2453 call alloc_compar_arrays(nfrag(1),nlevel)
2455 call card_concat(controlcard,.true.)
2456 write (iout,*) controlcard(:ilen(controlcard))
2457 call readi(controlcard,'NPIECE',npiece(j,1),0)
2458 call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
2459 call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
2460 call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
2461 call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
2462 call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
2463 call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
2464 call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
2465 call readi(controlcard,'RMS',irms(j,1),0)
2466 call readi(controlcard,'LOCAL',iloc(j),1)
2467 call readi(controlcard,'ELCONT',ielecont(j,1),1)
2468 if (ielecont(j,1).eq.0) then
2469 call readi(controlcard,'SCCONT',isccont(j,1),1)
2471 ang_cut(j)=ang_cut(j)*deg2rad
2472 ang_cut1(j)=ang_cut1(j)*deg2rad
2474 call card_concat(controlcard,.true.)
2475 call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
2476 call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
2478 write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
2479 (ifrag(1,k,j),ifrag(2,k,j),&
2480 k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
2481 " ang_cut1",ang_cut1(j)*rad2deg
2482 write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
2483 write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
2484 write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
2485 " ilocal",iloc(j)," isccont",isccont(j,1)
2487 ! Read data pertaning to higher levels
2489 call card_concat(controlcard,.true.)
2490 call readi(controlcard,'NFRAG',NFRAG(i),0)
2491 write (iout,*) "i",i," nfrag",nfrag(i)
2493 call card_concat(controlcard,.true.)
2495 call readi(controlcard,'ELCONT',ielecont(j,i),0)
2496 if (ielecont(j,i).eq.0) then
2497 call readi(controlcard,'SCCONT',isccont(j,i),1)
2499 call readi(controlcard,'RMS',irms(j,i),0)
2505 call readi(controlcard,'NPIECE',npiece(j,i),0)
2506 call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
2507 call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
2508 call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
2510 call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
2511 call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
2512 write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
2513 n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
2514 " isccont",isccont(j,i)," irms",irms(j,i)
2515 write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
2516 write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
2517 write(iout,*)"nc_frac",nc_fragm(j,i),&
2518 " nc_req",nc_req_setf(j,i)
2521 if (binary) write (iout,*) "Classes written in binary format."
2524 call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
2525 call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
2526 call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
2527 call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
2528 call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
2529 call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
2530 call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
2531 call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
2532 call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
2533 call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
2534 call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
2535 call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
2536 call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
2537 call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
2538 call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
2539 call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
2540 call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
2541 call readi(controlcard,'RMS_SINGLE',irms_single,0)
2542 call readi(controlcard,'CONT_SINGLE',icont_single,1)
2543 call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
2544 call readi(controlcard,'RMS_PAIR',irms_pair,0)
2545 call readi(controlcard,'CONT_PAIR',icont_pair,1)
2546 call readi(controlcard,'SPLIT_BET',isplit_bet,0)
2547 angcut_hel=angcut_hel*deg2rad
2548 angcut1_hel=angcut1_hel*deg2rad
2549 angcut_bet=angcut_bet*deg2rad
2550 angcut1_bet=angcut1_bet*deg2rad
2551 angcut_strand=angcut_strand*deg2rad
2552 angcut1_strand=angcut1_strand*deg2rad
2553 write (iout,*) "Automatic detection of structural elements"
2554 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2555 ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2556 ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
2557 ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2558 ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
2559 ' SPLIT_BET',isplit_bet
2560 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2561 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2562 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2563 ' MAXANG_HEL',angcut1_hel*rad2deg
2564 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2565 ' MAXANG_BET',angcut1_bet*rad2deg
2566 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2567 ' MAXANG_STRAND',angcut1_strand*rad2deg
2568 write (iout,*) 'FRAC_MIN',frac_min_set
2570 end subroutine read_compar
2571 !--------------------------------------------------------------------------------
2573 !--------------------------------------------------------------------------------
2574 subroutine read_ref_structure(*)
2576 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
2579 use control_data, only:pdbref
2580 use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
2581 nstart_sup,nstart_seq,nperm,nres0
2582 use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
2583 use compare, only:seq_comp !,contact,elecont
2584 use geometry, only:chainbuild,dist
2585 use io_config, only:readpdb
2587 use conform_compar, only:contact,elecont
2589 ! include 'DIMENSIONS'
2590 ! include 'DIMENSIONS.ZSCOPT'
2591 ! include 'DIMENSIONS.COMPAR'
2592 ! include 'COMMON.IOUNITS'
2593 ! include 'COMMON.GEO'
2594 ! include 'COMMON.VAR'
2595 ! include 'COMMON.INTERACT'
2596 ! include 'COMMON.LOCAL'
2597 ! include 'COMMON.NAMES'
2598 ! include 'COMMON.CHAIN'
2599 ! include 'COMMON.FFIELD'
2600 ! include 'COMMON.SBRIDGE'
2601 ! include 'COMMON.HEADER'
2602 ! include 'COMMON.CONTROL'
2603 ! include 'COMMON.CONTACTS1'
2604 ! include 'COMMON.PEPTCONT'
2605 ! include 'COMMON.TIME1'
2606 ! include 'COMMON.COMPAR'
2607 character(len=4) :: sequence(nres)
2609 !el real(kind=8) :: x(maxvar)
2610 integer :: itype_pdb(nres,5)
2611 !el logical seq_comp
2612 integer :: i,j,k,nres_pdb,iaux,mnum
2613 real(kind=8) :: ddsc !el,dist
2614 integer :: kkk !,ilen
2618 write (iout,*) "pdbref",pdbref
2620 read(inp,'(a)') pdbfile
2621 write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
2622 pdbfile(:ilen(pdbfile))
2623 open(ipdbin,file=pdbfile,status='old',err=33)
2625 33 write (iout,'(a)') 'Error opening PDB file.'
2630 itype_pdb(i,mnum)=itype(i,mnum)
2636 iaux=itype_pdb(i,mnum)
2637 itype_pdb(i,mnum)=itype(i,mnum)
2645 if (nsup.le.(nct-nnt+1)) then
2646 do i=0,nct-nnt+1-nsup
2647 if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
2649 do j=nnt+nsup-1,nnt,-1
2651 cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
2654 do j=nnt+nsup-1,nnt,-1
2656 cref(k,j+i,kkk)=cref(k,j,kkk)
2658 phi_ref(j+i)=phi_ref(j)
2659 theta_ref(j+i)=theta_ref(j)
2660 alph_ref(j+i)=alph_ref(j)
2661 omeg_ref(j+i)=omeg_ref(j)
2665 write (iout,'(i5,3f10.5,5x,3f10.5)') &
2666 j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
2674 write (iout,'(a)') &
2675 'Error - sequences to be superposed do not match.'
2678 do i=0,nsup-(nct-nnt+1)
2679 if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
2682 nstart_sup=nstart_sup+i
2687 write (iout,'(a)') &
2688 'Error - sequences to be superposed do not match.'
2692 write (iout,'(a,i5)') &
2693 'Experimental structure begins at residue',nstart_seq
2695 call read_angles(inp,*38)
2697 38 write (iout,'(a)') 'Error reading reference structure.'
2706 cref(j,i,kkk)=c(j,i)
2710 nend_sup=nstart_sup+nsup-1
2713 c(j,i)=cref(j,i,kkk)
2719 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
2721 if (itype(i,mnum).ne.10) then
2722 ddsc = dist(i,nres+i)
2724 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
2728 dc_norm(j,nres+i)=0.0d0
2731 ! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
2732 ! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
2733 ! dc_norm(3,nres+i)**2
2735 dc(j,i)=c(j,i+1)-c(j,i)
2739 dc_norm(j,i)=dc(j,i)/ddsc
2742 ! print *,"Calling contact"
2743 call contact(.true.,ncont_ref,icont_ref(1,1),&
2744 nstart_sup,nend_sup)
2745 ! print *,"Calling elecont"
2746 call elecont(.true.,ncont_pept_ref,&
2747 icont_pept_ref(1,1),&
2748 nstart_sup,nend_sup)
2749 write (iout,'(a,i3,a,i3,a,i3,a)') &
2750 'Number of residues to be superposed:',nsup,&
2751 ' (from residue',nstart_sup,' to residue',&
2754 end subroutine read_ref_structure
2755 !--------------------------------------------------------------------------------
2757 !--------------------------------------------------------------------------------
2758 subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
2760 use geometry_data, only:nres,c
2761 use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
2762 ! implicit real*8 (a-h,o-z)
2763 ! include 'DIMENSIONS'
2764 ! include 'DIMENSIONS.ZSCOPT'
2765 ! include 'COMMON.CHAIN'
2766 ! include 'COMMON.INTERACT'
2767 ! include 'COMMON.NAMES'
2768 ! include 'COMMON.IOUNITS'
2769 ! include 'COMMON.HEADER'
2770 ! include 'COMMON.SBRIDGE'
2771 character(len=50) :: tytul
2772 character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
2773 'D','E','F','G','H','I','J'/),shape(chainid))
2774 integer,dimension(nres) :: ica !(maxres)
2775 real(kind=8) :: temp,efree,etot,entropy,rmsdev
2776 integer :: ii,i,j,iti,ires,iatom,ichain,mnum
2777 write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
2779 write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
2781 write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
2789 if (iti.eq.ntyp1) then
2792 write (ipdb,'(a)') 'TER'
2797 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
2801 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
2802 ires,(c(j,nres+i),j=1,3)
2806 write (ipdb,'(a)') 'TER'
2809 if (itype(i,mnum).eq.ntyp1) cycle
2810 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
2811 write (ipdb,30) ica(i),ica(i+1)
2812 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
2813 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2814 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
2815 write (ipdb,30) ica(i),ica(i)+1
2818 if (itype(nct,molnum(nct)).ne.10) then
2819 write (ipdb,30) ica(nct),ica(nct)+1
2822 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2824 write (ipdb,'(a6)') 'ENDMDL'
2825 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2826 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2827 30 FORMAT ('CONECT',8I5)
2829 end subroutine pdboutW
2831 !------------------------------------------------------------------------------
2833 !-----------------------------------------------------------------------------
2834 !-----------------------------------------------------------------------------