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('BONDPAR_NUCL',bondname_nucl)
55 open (ibond_nucl,file=bondname_nucl,status='old')
56 call mygetenv('THETPAR',thetname)
57 open (ithep,file=thetname,status='old')
58 call mygetenv('ROTPAR',rotname)
59 open (irotam,file=rotname,status='old')
60 call mygetenv('TORPAR',torname)
61 open (itorp,file=torname,status='old')
62 call mygetenv('TORDPAR',tordname)
63 open (itordp,file=tordname,status='old')
64 call mygetenv('FOURIER',fouriername)
65 open (ifourier,file=fouriername,status='old')
66 call mygetenv('SCCORPAR',sccorname)
67 open (isccor,file=sccorname,status='old')
68 call mygetenv('ELEPAR',elename)
69 open (ielep,file=elename,status='old')
70 call mygetenv('SIDEPAR',sidename)
71 open (isidep,file=sidename,status='old')
72 call mygetenv('SIDEP',sidepname)
73 open (isidep1,file=sidepname,status="old")
74 call mygetenv('THETPAR_NUCL',thetname_nucl)
75 open (ithep_nucl,file=thetname_nucl,status='old')
76 call mygetenv('ROTPAR_NUCL',rotname_nucl)
77 open (irotam_nucl,file=rotname_nucl,status='old')
78 call mygetenv('TORPAR_NUCL',torname_nucl)
79 open (itorp_nucl,file=torname_nucl,status='old')
80 call mygetenv('TORDPAR_NUCL',tordname_nucl)
81 open (itordp_nucl,file=tordname_nucl,status='old')
82 call mygetenv('SIDEPAR_NUCL',sidename_nucl)
83 open (isidep_nucl,file=sidename_nucl,status='old')
84 call mygetenv('SIDEPAR_SCBASE',sidename_scbase)
85 open (isidep_scbase,file=sidename_scbase,status='old')
86 call mygetenv('PEPPAR_PEPBASE',pepname_pepbase)
87 open (isidep_pepbase,file=pepname_pepbase,status='old')
88 call mygetenv('SCPAR_PHOSPH',pepname_scpho)
89 open (isidep_scpho,file=pepname_scpho,status='old')
90 call mygetenv('PEPPAR_PHOSPH',pepname_peppho)
91 open (isidep_peppho,file=pepname_peppho,status='old')
92 call mygetenv('LIPTRANPAR',liptranname)
93 open (iliptranpar,file=liptranname,status='old',action='read')
94 call mygetenv('TUBEPAR',tubename)
95 open (itube,file=tubename,status='old',action='read')
96 call mygetenv('IONPAR',ionname)
97 open (iion,file=ionname,status='old',action='read')
99 call mygetenv('SCPPAR_NUCL',scpname_nucl)
100 open (iscpp_nucl,file=scpname_nucl,status='old')
105 ! 8/9/01 In the newest version SCp interaction constants are read from a file
106 ! Use -DOLDSCP to use hard-coded constants instead.
108 call mygetenv('SCPPAR',scpname)
109 open (iscpp,file=scpname,status='old')
112 if (MyID.eq.BossID) then
113 MyRank = MyID/fgProcs
116 print *,'OpenUnits: processor',MyRank
117 call numstr(MyRank,liczba)
118 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//liczba
120 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
122 open(iout,file=outname,status='unknown')
123 write (iout,'(80(1h-))')
124 write (iout,'(30x,a)') "FILE ASSIGNMENT"
125 write (iout,'(80(1h-))')
126 write (iout,*) "Input file : ",&
127 prefix(:ilen(prefix))//'.inp'
128 write (iout,*) "Output file : ",&
129 outname(:ilen(outname))
131 write (iout,*) "Sidechain potential file : ",&
132 sidename(:ilen(sidename))
134 write (iout,*) "SCp potential file : ",&
135 scpname(:ilen(scpname))
137 write (iout,*) "Electrostatic potential file : ",&
138 elename(:ilen(elename))
139 write (iout,*) "Cumulant coefficient file : ",&
140 fouriername(:ilen(fouriername))
141 write (iout,*) "Torsional parameter file : ",&
142 torname(:ilen(torname))
143 write (iout,*) "Double torsional parameter file : ",&
144 tordname(:ilen(tordname))
145 write (iout,*) "Backbone-rotamer parameter file : ",&
146 sccorname(:ilen(sccorname))
147 write (iout,*) "Bond & inertia constant file : ",&
148 bondname(:ilen(bondname))
149 write (iout,*) "Bending parameter file : ",&
150 thetname(:ilen(thetname))
151 write (iout,*) "Rotamer parameter file : ",&
152 rotname(:ilen(rotname))
153 write (iout,'(80(1h-))')
156 end subroutine openunits
157 !-----------------------------------------------------------------------------
159 !-----------------------------------------------------------------------------
160 subroutine molread(*)
162 ! Read molecular data.
165 use geometry_data, only:nres,deg2rad,c,dc,nres_molec
166 use control_data, only:iscode
167 use io_base, only:rescode
168 use control, only:setup_var,init_int_table
169 use geometry, only:alloc_geo_arrays
170 use energy, only:alloc_ener_arrays
171 ! implicit real*8 (a-h,o-z)
172 ! include 'DIMENSIONS'
173 ! include 'DIMENSIONS.ZSCOPT'
174 ! include 'COMMON.IOUNITS'
175 ! include 'COMMON.GEO'
176 ! include 'COMMON.VAR'
177 ! include 'COMMON.INTERACT'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.NAMES'
180 ! include 'COMMON.CHAIN'
181 ! include 'COMMON.FFIELD'
182 ! include 'COMMON.SBRIDGE'
183 ! include 'COMMON.TORCNSTR'
184 ! include 'COMMON.CONTROL'
185 character(len=4),dimension(:),allocatable :: sequence !(nres)
186 !el integer :: rescode
187 !el real(kind=8) :: x(maxvar)
188 character(len=320) :: controlcard !,ucase
189 integer,dimension(nres,5) :: itype_pdb !(maxres)
190 integer :: i,j,i1,i2,it1,it2,mnum
191 real(kind=8) :: scalscp
192 !el logical :: seq_comp
193 write(iout,*) "START MOLREAD"
197 print *,"nres_molec, initial",nres_molec(i)
200 call card_concat(controlcard,.true.)
201 call reada(controlcard,'SCAL14',scal14,0.4d0)
202 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
203 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
204 call reada(controlcard,'TEMP0',temp0,300.0d0) !el
205 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
206 r0_corr=cutoff_corr-delt_corr
207 call readi(controlcard,"NRES",nres_molec(1),0)
208 call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
209 call readi(controlcard,"NRES_CAT",nres_molec(5),0)
212 nres=nres_molec(i)+nres
214 ! allocate(sequence(nres+1))
215 !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
216 call alloc_geo_arrays
217 call alloc_ener_arrays
218 ! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
219 !----------------------------
220 allocate(c(3,2*nres+2))
221 allocate(dc(3,0:2*nres+2))
222 allocate(itype(nres+2,5))
223 allocate(itel(nres+2))
224 if (.not. allocated(molnum)) allocate(molnum(nres+2))
240 !--------------------------
242 allocate(sequence(nres+1))
243 iscode=index(controlcard,"ONE_LETTER")
245 write (iout,*) "Error: no residues in molecule"
248 if (nres.gt.maxres) then
249 write (iout,*) "Error: too many residues",nres,maxres
251 write(iout,*) 'nres=',nres
253 ! Read sequence of the protein
254 if (iscode.gt.0) then
255 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
257 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
259 ! Convert sequence to numeric code
263 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
265 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
268 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
270 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
273 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
276 write (iout,*) "Numeric code:"
277 write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres)
281 if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
283 if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
287 else if (iabs(itype(i+1,mnum)).ne.20) then
289 else if (iabs(itype(i,mnum)).ne.20) then
296 write (iout,*) "ITEL"
299 write (iout,*) i,itype(i,mnum),itel(i)
304 if (with_dihed_constr) then
306 read (inp,*) ndih_constr
307 if (ndih_constr.gt.0) then
309 write (iout,*) 'FTORS',ftors
310 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
312 'There are',ndih_constr,' constraints on phi angles.'
314 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
317 phi0(i)=deg2rad*phi0(i)
318 drange(i)=deg2rad*drange(i)
326 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
327 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
328 write(iout,*) 'NNT=',NNT,' NCT=',NCT
332 write (iout,'(/a,i3,a)') 'The chain contains',ns,&
333 ' disulfide-bridging cysteines.'
334 write (iout,'(20i4)') (iss(i),i=1,ns)
335 write (iout,'(/a/)') 'Pre-formed links are:'
342 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
343 restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',&
344 dhpb(i),ebr,forcon(i)
349 end subroutine molread
350 !-----------------------------------------------------------------------------
352 !-----------------------------------------------------------------------------
353 subroutine parmread(iparm,*)
358 ! Read the parameters of the probability distributions of the virtual-bond
359 ! valence angles and the side chains and energy parameters.
365 use control_data, only: maxterm,maxlor,maxterm_sccor,& !maxtor
366 maxtermd_1,maxtermd_2,tor_mode,scelemode !,maxthetyp,maxthetyp1
370 use io_config, only: printmat
371 use control, only: getenv_loc
378 ! implicit real*8 (a-h,o-z)
379 ! include 'DIMENSIONS'
380 ! include 'DIMENSIONS.ZSCOPT'
381 ! include 'DIMENSIONS.FREE'
382 ! include 'COMMON.IOUNITS'
383 ! include 'COMMON.CHAIN'
384 ! include 'COMMON.INTERACT'
385 ! include 'COMMON.GEO'
386 ! include 'COMMON.LOCAL'
387 ! include 'COMMON.TORSION'
388 ! include 'COMMON.FFIELD'
389 ! include 'COMMON.NAMES'
390 ! include 'COMMON.SBRIDGE'
391 ! include 'COMMON.WEIGHTS'
392 ! include 'COMMON.ENEPS'
393 ! include 'COMMON.SCCOR'
394 ! include 'COMMON.SCROT'
395 ! include 'COMMON.FREE'
396 character(len=1) :: t1,t2,t3
397 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
398 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
399 logical :: lprint,SPLIT_FOURIERTOR
400 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
401 character(len=800) :: controlcard
402 character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
403 tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
407 character(len=16) :: key
408 integer :: iparm,nkcctyp
409 !el real(kind=8) :: ip,mp
410 real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
411 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm, &
412 epspeptube,epssctube,sigmapeptube, &
415 real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
417 integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj
418 integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
424 ! Set LPRINT=.TRUE. for debugging
425 dwa16=2.0d0**(1.0d0/6.0d0)
428 ! Assign virtual-bond length
431 vblinv2=vblinv*vblinv
433 call card_concat(controlcard,.true.)
436 allocate(ww(max_eneW))
438 key = wname(i)(:ilen(wname(i)))
439 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
442 write (iout,*) "iparm",iparm," myparm",myparm
443 ! If reading not own parameters, skip assignment
445 if (iparm.eq.myparm .or. .not.separate_parset) then
448 ! Setup weights for UNRES
473 allocate(weights(n_ene))
488 weights(15)=0 !wstrain !
489 weights(16)=0 !wvdwpp !
491 weights(18)=0 !scal14 !
497 call card_concat(controlcard,.false.)
499 ! Return if not own parameters
501 if (iparm.ne.myparm .and. separate_parset) return
503 call reads(controlcard,"BONDPAR",bondname_t,bondname)
504 open (ibond,file=bondname_t,status='old')
506 call reads(controlcard,"THETPAR",thetname_t,thetname)
507 open (ithep,file=thetname_t,status='old')
509 call reads(controlcard,"ROTPAR",rotname_t,rotname)
510 open (irotam,file=rotname_t,status='old')
512 call reads(controlcard,"TORPAR",torname_t,torname)
513 open (itorp,file=torname_t,status='old')
515 call reads(controlcard,"TORDPAR",tordname_t,tordname)
516 open (itordp,file=tordname_t,status='old')
518 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
519 open (isccor,file=sccorname_t,status='old')
521 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
522 open (ifourier,file=fouriername_t,status='old')
524 call reads(controlcard,"ELEPAR",elename_t,elename)
525 open (ielep,file=elename_t,status='old')
527 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
528 open (isidep,file=sidename_t,status='old')
530 call reads(controlcard,"SCPPAR",scpname_t,scpname)
531 open (iscpp,file=scpname_t,status='old')
533 call getenv_loc('IONPAR',ionname)
534 open (iion,file=ionname,status='old')
536 write (iout,*) "Parameter set:",iparm
537 write (iout,*) "Energy-term weights:"
539 write (iout,'(a16,f10.5)') wname(i),ww(i)
541 write (iout,*) "Sidechain potential file : ",&
542 sidename_t(:ilen(sidename_t))
544 write (iout,*) "SCp potential file : ",&
545 scpname_t(:ilen(scpname_t))
547 write (iout,*) "Electrostatic potential file : ",&
548 elename_t(:ilen(elename_t))
549 write (iout,*) "Cumulant coefficient file : ",&
550 fouriername_t(:ilen(fouriername_t))
551 write (iout,*) "Torsional parameter file : ",&
552 torname_t(:ilen(torname_t))
553 write (iout,*) "Double torsional parameter file : ",&
554 tordname_t(:ilen(tordname_t))
555 write (iout,*) "Backbone-rotamer parameter file : ",&
556 sccorname(:ilen(sccorname))
557 write (iout,*) "Bond & inertia constant file : ",&
558 bondname_t(:ilen(bondname_t))
559 write (iout,*) "Bending parameter file : ",&
560 thetname_t(:ilen(thetname_t))
561 write (iout,*) "Rotamer parameter file : ",&
562 rotname_t(:ilen(rotname_t))
565 ! Read the virtual-bond parameters, masses, and moments of inertia
566 ! and Stokes' radii of the peptide group and side chains
568 allocate(dsc(ntyp1)) !(ntyp1)
569 allocate(dsc_inv(ntyp1)) !(ntyp1)
570 allocate(nbondterm(ntyp)) !(ntyp)
571 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
572 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
573 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
574 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
575 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
577 !el allocate(msc(ntyp+1)) !(ntyp+1)
578 !el allocate(isc(ntyp+1)) !(ntyp+1)
579 !el allocate(restok(ntyp+1)) !(ntyp+1)
580 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
583 read (ibond,*) vbldp0,akp
586 read (ibond,*) vbldsc0(1,i),aksc(1,i)
587 dsc(i) = vbldsc0(1,i)
591 dsc_inv(i)=1.0D0/dsc(i)
595 read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
597 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
599 dsc(i) = vbldsc0(1,i)
603 dsc_inv(i)=1.0D0/dsc(i)
608 write(iout,'(/a/)')"Force constants virtual bonds:"
609 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
611 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
613 write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
614 vbldsc0(1,i),aksc(1,i),abond0(1,i)
616 write (iout,'(13x,3f10.5)') &
617 vbldsc0(j,i),aksc(j,i),abond0(j,i)
621 if (.not. allocated(msc)) allocate(msc(ntyp1,5))
622 if (.not. allocated(restok)) allocate(restok(ntyp1,5))
625 read(iion,*) msc(i,5),restok(i,5)
626 print *,msc(i,5),restok(i,5)
630 read (iion,*) ncatprotparm
631 allocate(catprm(ncatprotparm,4))
633 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
636 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
639 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2)
640 ! dsc(i) = vbldsc0_nucl(1,i)
644 ! dsc_inv(i)=1.0D0/dsc(i)
649 !----------------------------------------------------
650 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
651 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
652 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
653 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
654 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
655 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
661 athet(j,i,ichir1,ichir2)=0.0D0
662 bthet(j,i,ichir1,ichir2)=0.0D0
676 !elwrite(iout,*) "parmread kontrol"
680 ! Read the parameters of the probability distribution/energy expression
681 ! of the virtual-bond valence angles theta
684 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
685 (bthet(j,i,1,1),j=1,2)
686 read (ithep,*) (polthet(j,i),j=0,3)
687 !elwrite(iout,*) "parmread kontrol in cryst_theta"
688 read (ithep,*) (gthet(j,i),j=1,3)
689 !elwrite(iout,*) "parmread kontrol in cryst_theta"
690 read (ithep,*) theta0(i),sig0(i),sigc0(i)
692 !elwrite(iout,*) "parmread kontrol in cryst_theta"
694 !elwrite(iout,*) "parmread kontrol in cryst_theta"
696 athet(1,i,1,-1)=athet(1,i,1,1)
697 athet(2,i,1,-1)=athet(2,i,1,1)
698 bthet(1,i,1,-1)=-bthet(1,i,1,1)
699 bthet(2,i,1,-1)=-bthet(2,i,1,1)
700 athet(1,i,-1,1)=-athet(1,i,1,1)
701 athet(2,i,-1,1)=-athet(2,i,1,1)
702 bthet(1,i,-1,1)=bthet(1,i,1,1)
703 bthet(2,i,-1,1)=bthet(2,i,1,1)
705 !elwrite(iout,*) "parmread kontrol in cryst_theta"
708 athet(1,i,-1,-1)=athet(1,-i,1,1)
709 athet(2,i,-1,-1)=-athet(2,-i,1,1)
710 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
711 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
712 athet(1,i,-1,1)=athet(1,-i,1,1)
713 athet(2,i,-1,1)=-athet(2,-i,1,1)
714 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
715 bthet(2,i,-1,1)=bthet(2,-i,1,1)
716 athet(1,i,1,-1)=-athet(1,-i,1,1)
717 athet(2,i,1,-1)=athet(2,-i,1,1)
718 bthet(1,i,1,-1)=bthet(1,-i,1,1)
719 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
724 polthet(j,i)=polthet(j,-i)
727 gthet(j,i)=gthet(j,-i)
730 !elwrite(iout,*) "parmread kontrol in cryst_theta"
732 !elwrite(iout,*) "parmread kontrol in cryst_theta"
735 ! & 'Parameters of the virtual-bond valence angles:'
736 ! write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
737 ! & ' ATHETA0 ',' A1 ',' A2 ',
740 ! write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
741 ! & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
743 ! write (iout,'(/a/9x,5a/79(1h-))')
744 ! & 'Parameters of the expression for sigma(theta_c):',
745 ! & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
746 ! & ' ALPH3 ',' SIGMA0C '
748 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
749 ! & (polthet(j,i),j=0,3),sigc0(i)
751 ! write (iout,'(/a/9x,5a/79(1h-))')
752 ! & 'Parameters of the second gaussian:',
753 ! & ' THETA0 ',' SIGMA0 ',' G1 ',
756 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
757 ! & sig0(i),(gthet(j,i),j=1,3)
760 'Parameters of the virtual-bond valence angles:'
761 write (iout,'(/a/9x,5a/79(1h-))') &
762 'Coefficients of expansion',&
763 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
764 ' b1*10^1 ',' b2*10^1 '
766 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
767 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
768 (10*bthet(j,i,1,1),j=1,2)
770 write (iout,'(/a/9x,5a/79(1h-))') &
771 'Parameters of the expression for sigma(theta_c):',&
772 ' alpha0 ',' alph1 ',' alph2 ',&
773 ' alhp3 ',' sigma0c '
775 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
776 (polthet(j,i),j=0,3),sigc0(i)
778 write (iout,'(/a/9x,5a/79(1h-))') &
779 'Parameters of the second gaussian:',&
780 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
783 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
784 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
789 ! Read the parameters of Utheta determined from ab initio surfaces
790 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
792 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
793 ! write (iout,*) "tu dochodze"
794 IF (tor_mode.eq.0) THEN
795 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
796 read (ithep,*) nthetyp,ntheterm,ntheterm2,&
797 ntheterm3,nsingle,ndouble
798 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
800 !----------------------------------------------------
801 ! allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
802 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
803 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
804 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
805 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
806 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
807 !(maxtheterm,-maxthetyp1:maxthetyp1,&
808 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
809 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
810 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
811 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
812 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
813 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
814 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
815 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
816 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
817 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
818 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
819 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
820 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
821 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
822 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
823 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
824 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
827 read (ithep,*) (ithetyp(i),i=1,ntyp1)
829 ithetyp(i)=-ithetyp(-i)
831 ! write (iout,*) "tu dochodze"
832 aa0thet(:,:,:,:)=0.0d0
833 aathet(:,:,:,:,:)=0.0d0
834 bbthet(:,:,:,:,:,:)=0.0d0
835 ccthet(:,:,:,:,:,:)=0.0d0
836 ddthet(:,:,:,:,:,:)=0.0d0
837 eethet(:,:,:,:,:,:)=0.0d0
838 ffthet(:,:,:,:,:,:,:)=0.0d0
839 ggthet(:,:,:,:,:,:,:)=0.0d0
843 do j=-nthetyp,nthetyp
844 do k=-nthetyp,nthetyp
845 read (ithep,'(6a)') res1
846 read (ithep,*) aa0thet(i,j,k,iblock)
847 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
849 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
850 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
851 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
852 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
855 (((ffthet(llll,lll,ll,i,j,k,iblock),&
856 ffthet(lll,llll,ll,i,j,k,iblock),&
857 ggthet(llll,lll,ll,i,j,k,iblock),&
858 ggthet(lll,llll,ll,i,j,k,iblock),&
859 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
864 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
865 ! coefficients of theta-and-gamma-dependent terms are zero.
870 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
871 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
873 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
874 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
877 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
879 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
882 ! Substitution for D aminoacids from symmetry.
885 do j=-nthetyp,nthetyp
886 do k=-nthetyp,nthetyp
887 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
889 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
893 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
894 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
895 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
896 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
902 ffthet(llll,lll,ll,i,j,k,iblock)= &
903 ffthet(llll,lll,ll,-i,-j,-k,iblock)
904 ffthet(lll,llll,ll,i,j,k,iblock)= &
905 ffthet(lll,llll,ll,-i,-j,-k,iblock)
906 ggthet(llll,lll,ll,i,j,k,iblock)= &
907 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
908 ggthet(lll,llll,ll,i,j,k,iblock)= &
909 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
919 ! Control printout of the coefficients of virtual-bond-angle potentials
923 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
927 write (iout,'(//4a)') &
928 'Type ',onelett(i),onelett(j),onelett(k)
929 write (iout,'(//a,10x,a)') " l","a[l]"
930 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
931 write (iout,'(i2,1pe15.5)') &
932 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
934 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
935 "b",l,"c",l,"d",l,"e",l
937 write (iout,'(i2,4(1pe15.5))') m,&
938 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
939 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
943 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
944 "f+",l,"f-",l,"g+",l,"g-",l
947 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
948 ffthet(n,m,l,i,j,k,iblock),&
949 ffthet(m,n,l,i,j,k,iblock),&
950 ggthet(n,m,l,i,j,k,iblock),&
951 ggthet(m,n,l,i,j,k,iblock)
962 !C here will be the apropriate recalibrating for D-aminoacid
963 read (ithep,*) nthetyp
964 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
965 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
966 do i=-nthetyp+1,nthetyp-1
967 read (ithep,*) nbend_kcc_Tb(i)
968 do j=0,nbend_kcc_Tb(i)
969 read (ithep,*) ijunk,v1bend_chyb(j,i)
974 "Parameters of the valence-only potentials"
975 do i=-nthetyp+1,nthetyp-1
976 write (iout,'(2a)') "Type ",toronelet(i)
977 do k=0,nbend_kcc_Tb(i)
978 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
985 !--------------- Reading theta parameters for nucleic acid-------
986 read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,&
987 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
988 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
989 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
990 allocate(aa0thet_nucl(nthetyp_nucl+1,&
991 nthetyp_nucl+1,nthetyp_nucl+1))
992 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
993 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
994 nthetyp_nucl+1,nthetyp_nucl+1))
995 !(maxtheterm,-maxthetyp1:maxthetyp1,&
996 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
997 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
998 nthetyp_nucl+1,nthetyp_nucl+1))
999 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1000 nthetyp_nucl+1,nthetyp_nucl+1))
1001 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1002 nthetyp_nucl+1,nthetyp_nucl+1))
1003 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1004 nthetyp_nucl+1,nthetyp_nucl+1))
1005 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1006 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1007 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1008 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1009 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1010 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1012 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1013 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1015 read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1017 aa0thet_nucl(:,:,:)=0.0d0
1018 aathet_nucl(:,:,:,:)=0.0d0
1019 bbthet_nucl(:,:,:,:,:)=0.0d0
1020 ccthet_nucl(:,:,:,:,:)=0.0d0
1021 ddthet_nucl(:,:,:,:,:)=0.0d0
1022 eethet_nucl(:,:,:,:,:)=0.0d0
1023 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1024 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1029 read (ithep_nucl,'(3a)') t1,t2,t3
1030 read (ithep_nucl,*) aa0thet_nucl(i,j,k)
1031 read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1032 read (ithep_nucl,*) &
1033 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1034 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1035 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1036 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1037 read (ithep_nucl,*) &
1038 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1039 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1040 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1046 !-------------------------------------------
1047 allocate(nlob(ntyp1)) !(ntyp1)
1048 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1049 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1050 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1068 gaussc(l,k,j,i)=0.0D0
1076 ! Read the parameters of the probability distribution/energy expression
1077 ! of the side chains.
1080 !c write (iout,*) "tu dochodze",i
1081 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
1085 dsc_inv(i)=1.0D0/dsc(i)
1096 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
1097 censc(1,1,-i)=censc(1,1,i)
1098 censc(2,1,-i)=censc(2,1,i)
1099 censc(3,1,-i)=-censc(3,1,i)
1101 read (irotam,*) bsc(j,i)
1102 read (irotam,*) (censc(k,j,i),k=1,3),&
1103 ((blower(k,l,j),l=1,k),k=1,3)
1104 censc(1,j,-i)=censc(1,j,i)
1105 censc(2,j,-i)=censc(2,j,i)
1106 censc(3,j,-i)=-censc(3,j,i)
1107 ! BSC is amplitude of Gaussian
1114 akl=akl+blower(k,m,j)*blower(l,m,j)
1118 if (((k.eq.3).and.(l.ne.3)) &
1119 .or.((l.eq.3).and.(k.ne.3))) then
1120 gaussc(k,l,j,-i)=-akl
1121 gaussc(l,k,j,-i)=-akl
1123 gaussc(k,l,j,-i)=akl
1124 gaussc(l,k,j,-i)=akl
1133 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1136 if (nlobi.gt.0) then
1137 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1138 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1139 ! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1140 ! write (iout,'(a,f10.4,4(16x,f10.4))')
1141 ! & 'Center ',(bsc(j,i),j=1,nlobi)
1142 ! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
1143 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1144 'log h',(bsc(j,i),j=1,nlobi)
1145 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1146 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1147 ! write (iout,'(a)')
1153 ! blower(k,l,j)=gaussc(ind,j,i)
1158 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1159 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1166 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1167 ! added by Urszula Kozlowska 07/11/2007
1169 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1177 read(irotam,*) sc_parmin(j,i)
1183 !---------reading nucleic acid parameters for rotamers-------------------
1184 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1185 do i=1,ntyp_molec(2)
1186 read (irotam_nucl,*)
1188 read(irotam_nucl,*) sc_parmin_nucl(j,i)
1194 write (iout,*) "Base rotamer parameters"
1195 do i=1,ntyp_molec(2)
1196 write (iout,'(a)') restyp(i,2)
1197 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1202 read (ifourier,*) nloctyp
1203 !el write(iout,*)"nloctyp",nloctyp
1204 SPLIT_FOURIERTOR = nloctyp.lt.0
1205 nloctyp = iabs(nloctyp)
1207 if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1))
1208 print *,"shape",shape(itype2loc)
1209 allocate(iloctyp(-nloctyp:nloctyp))
1210 allocate(bnew1(3,2,-nloctyp:nloctyp))
1211 allocate(bnew2(3,2,-nloctyp:nloctyp))
1212 allocate(ccnew(3,2,-nloctyp:nloctyp))
1213 allocate(ddnew(3,2,-nloctyp:nloctyp))
1214 allocate(e0new(3,-nloctyp:nloctyp))
1215 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1216 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1217 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1218 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1219 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1220 allocate(e0newtor(3,-nloctyp:nloctyp))
1221 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1223 read (ifourier,*) (itype2loc(i),i=1,ntyp)
1224 read (ifourier,*) (iloctyp(i),i=0,nloctyp-1)
1225 itype2loc(ntyp1)=nloctyp
1226 iloctyp(nloctyp)=ntyp1
1228 itype2loc(-i)=-itype2loc(i)
1231 allocate(iloctyp(-nloctyp:nloctyp))
1232 allocate(itype2loc(-ntyp1:ntyp1))
1239 iloctyp(-i)=-iloctyp(i)
1241 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1242 !c write (iout,*) "nloctyp",nloctyp,
1243 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1244 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1245 !c write (iout,*) "nloctyp",nloctyp,
1246 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1249 !c write (iout,*) "NEWCORR",i
1253 read (ifourier,*) bnew1(ii,j,i)
1256 !c write (iout,*) "NEWCORR BNEW1"
1257 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1260 read (ifourier,*) bnew2(ii,j,i)
1263 !c write (iout,*) "NEWCORR BNEW2"
1264 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1266 read (ifourier,*) ccnew(kk,1,i)
1267 read (ifourier,*) ccnew(kk,2,i)
1269 !c write (iout,*) "NEWCORR CCNEW"
1270 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1272 read (ifourier,*) ddnew(kk,1,i)
1273 read (ifourier,*) ddnew(kk,2,i)
1275 !c write (iout,*) "NEWCORR DDNEW"
1276 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1280 read (ifourier,*) eenew(ii,jj,kk,i)
1284 !c write (iout,*) "NEWCORR EENEW1"
1285 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1287 read (ifourier,*) e0new(ii,i)
1289 !c write (iout,*) (e0new(ii,i),ii=1,3)
1291 !c write (iout,*) "NEWCORR EENEW"
1294 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1295 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1296 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1297 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1302 bnew1(ii,1,-i)= bnew1(ii,1,i)
1303 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1304 bnew2(ii,1,-i)= bnew2(ii,1,i)
1305 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1308 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1309 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1310 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1311 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1312 ccnew(ii,1,-i)=ccnew(ii,1,i)
1313 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1314 ddnew(ii,1,-i)=ddnew(ii,1,i)
1315 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1317 e0new(1,-i)= e0new(1,i)
1318 e0new(2,-i)=-e0new(2,i)
1319 e0new(3,-i)=-e0new(3,i)
1321 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1322 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1323 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1324 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1328 write (iout,'(a)') "Coefficients of the multibody terms"
1329 do i=-nloctyp+1,nloctyp-1
1330 write (iout,*) "Type: ",onelet(iloctyp(i))
1331 write (iout,*) "Coefficients of the expansion of B1"
1333 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1335 write (iout,*) "Coefficients of the expansion of B2"
1337 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1339 write (iout,*) "Coefficients of the expansion of C"
1340 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1341 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1342 write (iout,*) "Coefficients of the expansion of D"
1343 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1344 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1345 write (iout,*) "Coefficients of the expansion of E"
1346 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1349 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1354 IF (SPLIT_FOURIERTOR) THEN
1356 !c write (iout,*) "NEWCORR TOR",i
1360 read (ifourier,*) bnew1tor(ii,j,i)
1363 !c write (iout,*) "NEWCORR BNEW1 TOR"
1364 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1367 read (ifourier,*) bnew2tor(ii,j,i)
1370 !c write (iout,*) "NEWCORR BNEW2 TOR"
1371 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1373 read (ifourier,*) ccnewtor(kk,1,i)
1374 read (ifourier,*) ccnewtor(kk,2,i)
1376 !c write (iout,*) "NEWCORR CCNEW TOR"
1377 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1379 read (ifourier,*) ddnewtor(kk,1,i)
1380 read (ifourier,*) ddnewtor(kk,2,i)
1382 !c write (iout,*) "NEWCORR DDNEW TOR"
1383 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1387 read (ifourier,*) eenewtor(ii,jj,kk,i)
1391 !c write (iout,*) "NEWCORR EENEW1 TOR"
1392 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1394 read (ifourier,*) e0newtor(ii,i)
1396 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1398 !c write (iout,*) "NEWCORR EENEW TOR"
1401 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1402 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1403 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1404 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1409 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1410 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1411 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1412 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1415 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1416 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1417 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1418 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1420 e0newtor(1,-i)= e0newtor(1,i)
1421 e0newtor(2,-i)=-e0newtor(2,i)
1422 e0newtor(3,-i)=-e0newtor(3,i)
1424 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1425 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1426 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1427 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1431 write (iout,'(a)') &
1432 "Single-body coefficients of the torsional potentials"
1433 do i=-nloctyp+1,nloctyp-1
1434 write (iout,*) "Type: ",onelet(iloctyp(i))
1435 write (iout,*) "Coefficients of the expansion of B1tor"
1437 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1438 j,(bnew1tor(k,j,i),k=1,3)
1440 write (iout,*) "Coefficients of the expansion of B2tor"
1442 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1443 j,(bnew2tor(k,j,i),k=1,3)
1445 write (iout,*) "Coefficients of the expansion of Ctor"
1446 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1447 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1448 write (iout,*) "Coefficients of the expansion of Dtor"
1449 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1450 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1451 write (iout,*) "Coefficients of the expansion of Etor"
1452 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1455 write (iout,'(1hE,2i1,2f10.5)') &
1456 j,k,(eenewtor(l,j,k,i),l=1,2)
1462 do i=-nloctyp+1,nloctyp-1
1465 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1470 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1474 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1475 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1476 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1477 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1480 ENDIF !SPLIT_FOURIER_TOR
1482 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1483 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1484 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1485 allocate(b(13,-nloctyp-1:nloctyp+1))
1487 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1490 read (ifourier,*) (b(ii,i),ii=1,13)
1492 write (iout,*) 'Type ',onelet(iloctyp(i))
1493 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1501 CCold(1,1,i)= b(7,i)
1502 CCold(2,2,i)=-b(7,i)
1503 CCold(2,1,i)= b(9,i)
1504 CCold(1,2,i)= b(9,i)
1505 CCold(1,1,-i)= b(7,i)
1506 CCold(2,2,-i)=-b(7,i)
1507 CCold(2,1,-i)=-b(9,i)
1508 CCold(1,2,-i)=-b(9,i)
1509 DDold(1,1,i)= b(6,i)
1510 DDold(2,2,i)=-b(6,i)
1511 DDold(2,1,i)= b(8,i)
1512 DDold(1,2,i)= b(8,i)
1513 DDold(1,1,-i)= b(6,i)
1514 DDold(2,2,-i)=-b(6,i)
1515 DDold(2,1,-i)=-b(8,i)
1516 DDold(1,2,-i)=-b(8,i)
1517 EEold(1,1,i)= b(10,i)+b(11,i)
1518 EEold(2,2,i)=-b(10,i)+b(11,i)
1519 EEold(2,1,i)= b(12,i)-b(13,i)
1520 EEold(1,2,i)= b(12,i)+b(13,i)
1521 EEold(1,1,-i)= b(10,i)+b(11,i)
1522 EEold(2,2,-i)=-b(10,i)+b(11,i)
1523 EEold(2,1,-i)=-b(12,i)+b(13,i)
1524 EEold(1,2,-i)=-b(12,i)-b(13,i)
1525 write(iout,*) "TU DOCHODZE"
1531 "Coefficients of the cumulants (independent of valence angles)"
1532 do i=-nloctyp+1,nloctyp-1
1533 write (iout,*) 'Type ',onelet(iloctyp(i))
1535 write(iout,'(2f10.5)') B(3,i),B(5,i)
1537 write(iout,'(2f10.5)') B(2,i),B(4,i)
1540 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1544 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1548 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1555 ! Read torsional parameters in old format
1557 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1559 read (itorp,*) ntortyp,nterm_old
1560 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1561 read (itorp,*) (itortyp(i),i=1,ntyp)
1563 !el from energy module--------
1564 allocate(v1(nterm_old,ntortyp,ntortyp))
1565 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1566 !el---------------------------
1572 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1578 write (iout,'(/a/)') 'Torsional constants:'
1581 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1582 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1590 ! Read torsional parameters
1592 IF (TOR_MODE.eq.0) THEN
1594 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1596 read (itorp,*) ntortyp
1597 read (itorp,*) (itortyp(i),i=1,ntyp)
1598 write (iout,*) 'ntortyp',ntortyp
1600 !el from energy module---------
1601 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1602 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1604 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1605 allocate(vlor2(maxlor,ntortyp,ntortyp))
1606 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1607 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1609 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1610 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1611 !el---------------------------
1613 do i=-ntortyp,ntortyp
1614 do j=-ntortyp,ntortyp
1620 !el---------------------------
1624 itortyp(i)=-itortyp(-i)
1626 ! write (iout,*) 'ntortyp',ntortyp
1628 do j=-ntortyp+1,ntortyp-1
1629 read (itorp,*) nterm(i,j,iblock),&
1631 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1632 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1635 do k=1,nterm(i,j,iblock)
1636 read (itorp,*) kk,v1(k,i,j,iblock),&
1638 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1639 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1640 v0ij=v0ij+si*v1(k,i,j,iblock)
1643 do k=1,nlor(i,j,iblock)
1644 read (itorp,*) kk,vlor1(k,i,j),&
1645 vlor2(k,i,j),vlor3(k,i,j)
1646 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1649 v0(-i,-j,iblock)=v0ij
1656 write (iout,'(/a/)') 'Torsional constants:'
1659 write (iout,*) 'ityp',i,' jtyp',j
1660 write (iout,*) 'Fourier constants'
1661 do k=1,nterm(i,j,iblock)
1662 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1665 write (iout,*) 'Lorenz constants'
1666 do k=1,nlor(i,j,iblock)
1667 write (iout,'(3(1pe15.5))') &
1668 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1675 ! 6/23/01 Read parameters for double torsionals
1677 !el from energy module------------
1678 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1679 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1680 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1681 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1682 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1683 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1684 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1685 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1686 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1687 !---------------------------------
1691 do j=-ntortyp+1,ntortyp-1
1692 do k=-ntortyp+1,ntortyp-1
1693 read (itordp,'(3a1)') t1,t2,t3
1694 ! write (iout,*) "OK onelett",
1697 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1698 .or. t3.ne.toronelet(k)) then
1699 write (iout,*) "Error in double torsional parameter file",&
1702 call MPI_Finalize(Ierror)
1704 stop "Error in double torsional parameter file"
1706 read (itordp,*) ntermd_1(i,j,k,iblock),&
1707 ntermd_2(i,j,k,iblock)
1708 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1709 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1710 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1711 ntermd_1(i,j,k,iblock))
1712 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1713 ntermd_1(i,j,k,iblock))
1714 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1715 ntermd_1(i,j,k,iblock))
1716 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1717 ntermd_1(i,j,k,iblock))
1718 ! Martix of D parameters for one dimesional foureir series
1719 do l=1,ntermd_1(i,j,k,iblock)
1720 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1721 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1722 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1723 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1724 ! write(iout,*) "whcodze" ,
1725 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1727 read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1728 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1729 v2s(m,l,i,j,k,iblock),&
1730 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1731 ! Martix of D parameters for two dimesional fourier series
1732 do l=1,ntermd_2(i,j,k,iblock)
1734 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1735 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1736 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1737 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1746 write (iout,*) 'Constants for double torsionals'
1749 do j=-ntortyp+1,ntortyp-1
1750 do k=-ntortyp+1,ntortyp-1
1751 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1752 ' nsingle',ntermd_1(i,j,k,iblock),&
1753 ' ndouble',ntermd_2(i,j,k,iblock)
1755 write (iout,*) 'Single angles:'
1756 do l=1,ntermd_1(i,j,k,iblock)
1757 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1758 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1759 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1760 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1763 write (iout,*) 'Pairs of angles:'
1764 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1765 do l=1,ntermd_2(i,j,k,iblock)
1766 write (iout,'(i5,20f10.5)') &
1767 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1770 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1771 do l=1,ntermd_2(i,j,k,iblock)
1772 write (iout,'(i5,20f10.5)') &
1773 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1774 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1784 itype2loc(i)=itortyp(i)
1787 ELSE IF (TOR_MODE.eq.1) THEN
1789 !C read valence-torsional parameters
1790 read (itorp,*) ntortyp
1792 write (iout,*) "Valence-torsional parameters read in ntortyp",&
1794 read (itorp,*) (itortyp(i),i=1,ntyp)
1795 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1798 itype2loc(i)=itortyp(i)
1802 itortyp(i)=-itortyp(-i)
1804 do i=-ntortyp+1,ntortyp-1
1805 do j=-ntortyp+1,ntortyp-1
1806 !C first we read the cos and sin gamma parameters
1807 read (itorp,'(13x,a)') string
1808 write (iout,*) i,j,string
1810 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1811 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1812 do k=1,nterm_kcc(j,i)
1813 do l=1,nterm_kcc_Tb(j,i)
1814 do ll=1,nterm_kcc_Tb(j,i)
1815 read (itorp,*) ii,jj,kk, &
1816 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1824 !c AL 4/8/16: Calculate coefficients from one-body parameters
1826 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1827 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
1828 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
1829 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1830 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1833 print *,i,itortyp(i)
1834 itortyp(i)=itype2loc(i)
1837 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
1839 do i=-ntortyp+1,ntortyp-1
1840 do j=-ntortyp+1,ntortyp-1
1843 do k=1,nterm_kcc_Tb(j,i)
1844 do l=1,nterm_kcc_Tb(j,i)
1845 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
1846 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1847 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
1848 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1851 do k=1,nterm_kcc_Tb(j,i)
1852 do l=1,nterm_kcc_Tb(j,i)
1854 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1855 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1856 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1857 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1859 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1860 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1861 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1862 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1866 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1870 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1874 if (tor_mode.gt.0 .and. lprint) then
1875 !c Print valence-torsional parameters
1876 write (iout,'(a)') &
1877 "Parameters of the valence-torsional potentials"
1878 do i=-ntortyp+1,ntortyp-1
1879 do j=-ntortyp+1,ntortyp-1
1880 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1881 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1882 do k=1,nterm_kcc(j,i)
1883 do l=1,nterm_kcc_Tb(j,i)
1884 do ll=1,nterm_kcc_Tb(j,i)
1885 write (iout,'(3i5,2f15.4)')&
1886 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1895 !elwrite(iout,*) "parmread kontrol sc-bb"
1896 ! Read of Side-chain backbone correlation parameters
1897 ! Modified 11 May 2012 by Adasko
1900 read (isccor,*) nsccortyp
1903 !c maxinter is maximum interaction sites
1904 !write(iout,*)"maxterm_sccor",maxterm_sccor
1905 !el from module energy-------------
1906 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1907 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1908 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1909 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1910 !-----------------------------------
1911 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1912 !-----------------------------------
1913 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1914 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1915 -nsccortyp:nsccortyp))
1916 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1917 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1918 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1919 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1920 !-----------------------------------
1921 do i=-nsccortyp,nsccortyp
1922 do j=-nsccortyp,nsccortyp
1926 !-----------------------------------
1928 read (isccor,*) (isccortyp(i),i=1,ntyp)
1930 isccortyp(i)=-isccortyp(-i)
1932 iscprol=isccortyp(20)
1933 ! write (iout,*) 'ntortyp',ntortyp
1935 !c maxinter is maximum interaction sites
1940 nterm_sccor(i,j),nlor_sccor(i,j)
1946 nterm_sccor(-i,j)=nterm_sccor(i,j)
1947 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1948 nterm_sccor(i,-j)=nterm_sccor(i,j)
1949 do k=1,nterm_sccor(i,j)
1950 read (isccor,*) kk,v1sccor(k,l,i,j),&
1952 if (j.eq.iscprol) then
1953 if (i.eq.isccortyp(10)) then
1954 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1955 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1957 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1958 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1959 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1960 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1961 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1962 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1963 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1964 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1967 if (i.eq.isccortyp(10)) then
1968 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1969 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1971 if (j.eq.isccortyp(10)) then
1972 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1973 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1975 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1976 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1977 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1978 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1979 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1980 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1984 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1985 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1986 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1987 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1990 do k=1,nlor_sccor(i,j)
1991 read (isccor,*) kk,vlor1sccor(k,i,j),&
1992 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1993 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1994 (1+vlor3sccor(k,i,j)**2)
1996 v0sccor(l,i,j)=v0ijsccor
1997 v0sccor(l,-i,j)=v0ijsccor1
1998 v0sccor(l,i,-j)=v0ijsccor2
1999 v0sccor(l,-i,-j)=v0ijsccor3
2005 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
2008 write (iout,*) 'ityp',i,' jtyp',j
2009 write (iout,*) 'Fourier constants'
2010 do k=1,nterm_sccor(i,j)
2011 write (iout,'(2(1pe15.5))') &
2012 (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
2014 write (iout,*) 'Lorenz constants'
2015 do k=1,nlor_sccor(i,j)
2016 write (iout,'(3(1pe15.5))') &
2017 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2024 ! Read electrostatic-interaction parameters
2027 write (iout,'(/a)') 'Electrostatic interaction constants:'
2028 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2029 'IT','JT','APP','BPP','AEL6','AEL3'
2031 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
2032 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
2033 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
2034 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
2039 app (i,j)=epp(i,j)*rri*rri
2040 bpp (i,j)=-2.0D0*epp(i,j)*rri
2041 ael6(i,j)=elpp6(i,j)*4.2D0**6
2042 ael3(i,j)=elpp3(i,j)*4.2D0**3
2043 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2048 ! Read side-chain interaction parameters.
2050 !el from module energy - COMMON.INTERACT-------
2051 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2052 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2053 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2054 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2055 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2056 allocate(epslip(ntyp,ntyp))
2067 !--------------------------------
2069 read (isidep,*) ipot,expon
2070 !el if (ipot.lt.1 .or. ipot.gt.5) then
2071 ! write (iout,'(2a)') 'Error while reading SC interaction',&
2072 ! 'potential file - unknown potential type.'
2076 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2077 ', exponents are ',expon,2*expon
2078 ! goto (10,20,30,30,40) ipot
2080 !----------------------- LJ potential ---------------------------------
2082 ! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2083 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2085 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2086 write (iout,'(a/)') 'The epsilon array:'
2087 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2088 write (iout,'(/a)') 'One-body parameters:'
2089 write (iout,'(a,4x,a)') 'residue','sigma'
2090 write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
2093 !----------------------- LJK potential --------------------------------
2095 ! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2096 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2097 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2099 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2100 write (iout,'(a/)') 'The epsilon array:'
2101 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2102 write (iout,'(/a)') 'One-body parameters:'
2103 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2104 write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2108 !---------------------- GB or BP potential -----------------------------
2111 if (scelemode.eq.0) then
2113 read (isidep,*)(eps(i,j),j=i,ntyp)
2115 read (isidep,*)(sigma0(i),i=1,ntyp)
2116 read (isidep,*)(sigii(i),i=1,ntyp)
2117 read (isidep,*)(chip(i),i=1,ntyp)
2118 read (isidep,*)(alp(i),i=1,ntyp)
2120 read (isidep,*)(epslip(i,j),j=i,ntyp)
2122 ! For the GB potential convert sigma'**2 into chi'
2125 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2129 write (iout,'(/a/)') 'Parameters of the BP potential:'
2130 write (iout,'(a/)') 'The epsilon array:'
2131 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2132 write (iout,'(/a)') 'One-body parameters:'
2133 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2135 write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
2136 chip(i),alp(i),i=1,ntyp)
2139 allocate(icharge(ntyp1))
2140 ! print *,ntyp,icharge(i)
2142 read (isidep,*) (icharge(i),i=1,ntyp)
2143 print *,ntyp,icharge(i)
2144 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2145 write (2,*) "icharge",(icharge(i),i=1,ntyp)
2146 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2147 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2148 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2149 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2150 allocate(epsintab(ntyp,ntyp))
2151 allocate(dtail(2,ntyp,ntyp))
2152 print *,"control line 1"
2153 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2154 allocate(wqdip(2,ntyp,ntyp))
2155 allocate(wstate(4,ntyp,ntyp))
2156 allocate(dhead(2,2,ntyp,ntyp))
2157 allocate(nstate(ntyp,ntyp))
2158 allocate(debaykap(ntyp,ntyp))
2159 print *,"control line 2"
2160 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2161 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2165 ! write (*,*) "Im in ALAB", i, " ", j
2167 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2168 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2169 chis(i,j),chis(j,i), &
2170 nstate(i,j),(wstate(k,i,j),k=1,4), &
2171 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2172 dtail(1,i,j),dtail(2,i,j), &
2173 epshead(i,j),sig0head(i,j), &
2174 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2175 alphapol(i,j),alphapol(j,i), &
2176 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
2177 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2183 sigma(i,j) = sigma(j,i)
2184 sigmap1(i,j)=sigmap1(j,i)
2185 sigmap2(i,j)=sigmap2(j,i)
2186 sigiso1(i,j)=sigiso1(j,i)
2187 sigiso2(i,j)=sigiso2(j,i)
2188 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2189 nstate(i,j) = nstate(j,i)
2190 dtail(1,i,j) = dtail(1,j,i)
2191 dtail(2,i,j) = dtail(2,j,i)
2193 alphasur(k,i,j) = alphasur(k,j,i)
2194 wstate(k,i,j) = wstate(k,j,i)
2195 alphiso(k,i,j) = alphiso(k,j,i)
2198 dhead(2,1,i,j) = dhead(1,1,j,i)
2199 dhead(2,2,i,j) = dhead(1,2,j,i)
2200 dhead(1,1,i,j) = dhead(2,1,j,i)
2201 dhead(1,2,i,j) = dhead(2,2,j,i)
2203 epshead(i,j) = epshead(j,i)
2204 sig0head(i,j) = sig0head(j,i)
2207 wqdip(k,i,j) = wqdip(k,j,i)
2210 wquad(i,j) = wquad(j,i)
2211 epsintab(i,j) = epsintab(j,i)
2212 debaykap(i,j)=debaykap(j,i)
2213 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2219 !--------------------- GBV potential -----------------------------------
2221 ! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2222 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2223 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2224 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2226 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2227 write (iout,'(a/)') 'The epsilon array:'
2228 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2229 write (iout,'(/a)') 'One-body parameters:'
2230 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2231 's||/s_|_^2',' chip ',' alph '
2232 write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2233 sigii(i),chip(i),alp(i),i=1,ntyp)
2236 write (iout,'(2a)') 'Error while reading SC interaction',&
2237 'potential file - unknown potential type.'
2243 !-----------------------------------------------------------------------
2244 ! Calculate the "working" parameters of SC interactions.
2246 !el from module energy - COMMON.INTERACT-------
2247 ! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2248 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1))
2249 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) !(ntyp,ntyp)
2250 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2251 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2252 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2253 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2255 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2263 if (scelemode.eq.0) then
2270 !--------------------------------
2275 epslip(i,j)=epslip(j,i)
2278 if (scelemode.eq.0) then
2281 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2282 sigma(j,i)=sigma(i,j)
2283 rs0(i,j)=dwa16*sigma(i,j)
2288 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2289 'Working parameters of the SC interactions:',&
2290 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2295 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2305 sigeps=dsign(1.0D0,epsij)
2307 aa_aq(i,j)=epsij*rrij*rrij
2308 bb_aq(i,j)=-sigeps*epsij*rrij
2309 aa_aq(j,i)=aa_aq(i,j)
2310 bb_aq(j,i)=bb_aq(i,j)
2311 epsijlip=epslip(i,j)
2312 sigeps=dsign(1.0D0,epsijlip)
2313 epsijlip=dabs(epsijlip)
2314 aa_lip(i,j)=epsijlip*rrij*rrij
2315 bb_lip(i,j)=-sigeps*epsijlip*rrij
2316 aa_lip(j,i)=aa_lip(i,j)
2317 bb_lip(j,i)=bb_lip(i,j)
2318 if ((ipot.gt.2).and. (scelemode.eq.0))then
2319 sigt1sq=sigma0(i)**2
2320 sigt2sq=sigma0(j)**2
2323 ratsig1=sigt2sq/sigt1sq
2324 ratsig2=1.0D0/ratsig1
2325 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2326 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2327 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2331 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2332 sigmaii(i,j)=rsum_max
2333 sigmaii(j,i)=rsum_max
2335 ! sigmaii(i,j)=r0(i,j)
2336 ! sigmaii(j,i)=r0(i,j)
2338 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2339 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2340 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2341 augm(i,j)=epsij*r_augm**(2*expon)
2342 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2349 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2350 restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2351 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2355 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2356 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2357 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2358 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2359 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2360 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2361 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2362 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2363 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2364 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2365 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2366 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2367 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2368 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2369 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2370 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2379 read (isidep_nucl,*) ipot_nucl
2380 ! print *,"TU?!",ipot_nucl
2381 if (ipot_nucl.eq.1) then
2382 do i=1,ntyp_molec(2)
2383 do j=i,ntyp_molec(2)
2384 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2385 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2389 do i=1,ntyp_molec(2)
2390 do j=i,ntyp_molec(2)
2391 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2392 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2393 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2397 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2398 do i=1,ntyp_molec(2)
2399 do j=i,ntyp_molec(2)
2400 rrij=sigma_nucl(i,j)
2404 epsij=4*eps_nucl(i,j)
2405 sigeps=dsign(1.0D0,epsij)
2407 aa_nucl(i,j)=epsij*rrij*rrij
2408 bb_nucl(i,j)=-sigeps*epsij*rrij
2409 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2410 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2411 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2412 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2413 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2414 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2417 aa_nucl(i,j)=aa_nucl(j,i)
2418 bb_nucl(i,j)=bb_nucl(j,i)
2419 ael3_nucl(i,j)=ael3_nucl(j,i)
2420 ael6_nucl(i,j)=ael6_nucl(j,i)
2421 ael63_nucl(i,j)=ael63_nucl(j,i)
2422 ael32_nucl(i,j)=ael32_nucl(j,i)
2423 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2424 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2425 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2426 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2427 eps_nucl(i,j)=eps_nucl(j,i)
2428 sigma_nucl(i,j)=sigma_nucl(j,i)
2429 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2433 write(iout,*) "tube param"
2434 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2435 ccavtubpep,dcavtubpep,tubetranenepep
2436 sigmapeptube=sigmapeptube**6
2437 sigeps=dsign(1.0D0,epspeptube)
2438 epspeptube=dabs(epspeptube)
2439 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2440 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2441 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2443 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2444 ccavtub(i),dcavtub(i),tubetranene(i)
2445 sigmasctube=sigmasctube**6
2446 sigeps=dsign(1.0D0,epssctube)
2447 epssctube=dabs(epssctube)
2448 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2449 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2450 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2452 !-----------------READING SC BASE POTENTIALS-----------------------------
2453 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2454 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2455 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2456 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2457 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2458 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2459 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2460 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2461 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2462 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2463 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2464 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2465 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2466 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2467 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2468 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2471 do i=1,ntyp_molec(1)
2472 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2473 write (*,*) "Im in ", i, " ", j
2474 read(isidep_scbase,*) &
2475 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2476 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2477 write(*,*) "eps",eps_scbase(i,j)
2478 read(isidep_scbase,*) &
2479 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2480 chis_scbase(i,j,1),chis_scbase(i,j,2)
2481 read(isidep_scbase,*) &
2482 dhead_scbasei(i,j), &
2483 dhead_scbasej(i,j), &
2484 rborn_scbasei(i,j),rborn_scbasej(i,j)
2485 read(isidep_scbase,*) &
2486 (wdipdip_scbase(k,i,j),k=1,3), &
2487 (wqdip_scbase(k,i,j),k=1,2)
2488 read(isidep_scbase,*) &
2489 alphapol_scbase(i,j), &
2490 epsintab_scbase(i,j)
2493 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2494 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2496 do i=1,ntyp_molec(1)
2497 do j=1,ntyp_molec(2)-1
2498 epsij=eps_scbase(i,j)
2499 rrij=sigma_scbase(i,j)
2504 sigeps=dsign(1.0D0,epsij)
2506 aa_scbase(i,j)=epsij*rrij*rrij
2507 bb_scbase(i,j)=-sigeps*epsij*rrij
2510 !-----------------READING PEP BASE POTENTIALS-------------------
2511 allocate(eps_pepbase(ntyp_molec(2)))
2512 allocate(sigma_pepbase(ntyp_molec(2)))
2513 allocate(chi_pepbase(ntyp_molec(2),2))
2514 allocate(chipp_pepbase(ntyp_molec(2),2))
2515 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2516 allocate(sigmap1_pepbase(ntyp_molec(2)))
2517 allocate(sigmap2_pepbase(ntyp_molec(2)))
2518 allocate(chis_pepbase(ntyp_molec(2),2))
2519 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2522 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2523 write (*,*) "Im in ", i, " ", j
2524 read(isidep_pepbase,*) &
2525 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2526 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2527 write(*,*) "eps",eps_pepbase(j)
2528 read(isidep_pepbase,*) &
2529 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2530 chis_pepbase(j,1),chis_pepbase(j,2)
2531 read(isidep_pepbase,*) &
2532 (wdipdip_pepbase(k,j),k=1,3)
2534 allocate(aa_pepbase(ntyp_molec(2)))
2535 allocate(bb_pepbase(ntyp_molec(2)))
2537 do j=1,ntyp_molec(2)-1
2538 epsij=eps_pepbase(j)
2539 rrij=sigma_pepbase(j)
2544 sigeps=dsign(1.0D0,epsij)
2546 aa_pepbase(j)=epsij*rrij*rrij
2547 bb_pepbase(j)=-sigeps*epsij*rrij
2549 !--------------READING SC PHOSPHATE-------------------------------------
2550 !--------------READING SC PHOSPHATE-------------------------------------
2551 allocate(eps_scpho(ntyp_molec(1)))
2552 allocate(sigma_scpho(ntyp_molec(1)))
2553 allocate(chi_scpho(ntyp_molec(1),2))
2554 allocate(chipp_scpho(ntyp_molec(1),2))
2555 allocate(alphasur_scpho(4,ntyp_molec(1)))
2556 allocate(sigmap1_scpho(ntyp_molec(1)))
2557 allocate(sigmap2_scpho(ntyp_molec(1)))
2558 allocate(chis_scpho(ntyp_molec(1),2))
2559 allocate(wqq_scpho(ntyp_molec(1)))
2560 allocate(wqdip_scpho(2,ntyp_molec(1)))
2561 allocate(alphapol_scpho(ntyp_molec(1)))
2562 allocate(epsintab_scpho(ntyp_molec(1)))
2563 allocate(dhead_scphoi(ntyp_molec(1)))
2564 allocate(rborn_scphoi(ntyp_molec(1)))
2565 allocate(rborn_scphoj(ntyp_molec(1)))
2566 allocate(alphi_scpho(ntyp_molec(1)))
2570 do j=1,ntyp_molec(1) ! without U then we will take T for U
2571 write (*,*) "Im in scpho ", i, " ", j
2572 read(isidep_scpho,*) &
2573 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2574 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2575 write(*,*) "eps",eps_scpho(j)
2576 read(isidep_scpho,*) &
2577 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2578 chis_scpho(j,1),chis_scpho(j,2)
2579 read(isidep_scpho,*) &
2580 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2581 read(isidep_scpho,*) &
2582 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2586 allocate(aa_scpho(ntyp_molec(1)))
2587 allocate(bb_scpho(ntyp_molec(1)))
2588 do j=1,ntyp_molec(1)
2595 sigeps=dsign(1.0D0,epsij)
2597 aa_scpho(j)=epsij*rrij*rrij
2598 bb_scpho(j)=-sigeps*epsij*rrij
2602 read(isidep_peppho,*) &
2603 eps_peppho,sigma_peppho
2604 read(isidep_peppho,*) &
2605 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2606 read(isidep_peppho,*) &
2607 (wqdip_peppho(k),k=1,2)
2615 sigeps=dsign(1.0D0,epsij)
2617 aa_peppho=epsij*rrij*rrij
2618 bb_peppho=-sigeps*epsij*rrij
2621 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2629 ! Define the SC-p interaction constants
2638 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2639 read (itorp_nucl,*) ntortyp_nucl
2640 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2641 !el from energy module---------
2642 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2643 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2645 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2646 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2647 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2648 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2650 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2651 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2652 !el---------------------------
2655 !el--------------------
2656 read (itorp_nucl,*) &
2657 (itortyp_nucl(i),i=1,ntyp_molec(2))
2658 ! print *,itortyp_nucl(:)
2659 !c write (iout,*) 'ntortyp',ntortyp
2662 read (itorp_nucl,*) nterm_nucl(i,j),nlor_nucl(i,j)
2663 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2666 do k=1,nterm_nucl(i,j)
2667 read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2668 v0ij=v0ij+si*v1_nucl(k,i,j)
2671 do k=1,nlor_nucl(i,j)
2672 read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),&
2673 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2674 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2681 !elwrite(iout,*) "parmread kontrol before oldscp"
2683 ! Define the SC-p interaction constants
2687 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2689 ! aad(i,1)=0.3D0*4.0D0**12
2690 ! Following line for constants currently implemented
2691 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2692 aad(i,1)=1.5D0*4.0D0**12
2693 ! aad(i,1)=0.17D0*5.6D0**12
2695 ! "Soft" SC-p repulsion
2697 ! Following line for constants currently implemented
2698 ! aad(i,1)=0.3D0*4.0D0**6
2699 ! "Hard" SC-p repulsion
2700 bad(i,1)=3.0D0*4.0D0**6
2701 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2710 ! 8/9/01 Read the SC-p interaction constants from file
2713 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
2716 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2717 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2718 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2719 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2723 write (iout,*) "Parameters of SC-p interactions:"
2725 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2726 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2730 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2732 do i=1,ntyp_molec(2)
2733 read (iscpp_nucl,*) eps_scp_nucl(i),rscp_nucl(i)
2735 do i=1,ntyp_molec(2)
2736 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2737 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2739 r0pp=1.12246204830937298142*5.16158
2745 ! Define the constants of the disulfide bridge
2749 ! Old arbitrary potential - commented out.
2754 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2755 ! energy surface of diethyl disulfide.
2756 ! A. Liwo and U. Kozlowska, 11/24/03
2767 write (iout,'(/a)') "Disulfide bridge parameters:"
2768 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2769 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2770 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2771 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2775 end subroutine parmread
2777 !-----------------------------------------------------------------------------
2779 !-----------------------------------------------------------------------------
2780 subroutine mygetenv(string,var)
2784 ! This subroutine passes the environmental variables to FORTRAN program.
2785 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
2786 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
2787 ! reads the environmental variables from $HOME/.env
2789 ! Usage: As for the standard FORTRAN GETENV subroutine.
2791 ! Purpose: some versions/installations of MPI do not transfer the environmental
2792 ! variables to slave processors, if these variables are set in the shell script
2793 ! from which mpirun is called.
2802 character*(*) :: string,var
2803 #if defined(MYGETENV) && defined(MPI)
2804 ! include "DIMENSIONS.ZSCOPT"
2806 ! include "COMMON.MPI"
2807 !el character*360 ucase
2809 character(len=360) :: string1(360),karta
2810 character(len=240) :: home
2813 call getenv("HOME",home)
2814 open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
2816 read (99,end=111,err=111,'(a)') karta
2820 call split_string(karta,string1,80,n)
2821 if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
2822 string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
2824 print *,"Processor",me,": ",var(:ilen(var)),&
2825 " assigned to ",string(:ilen(string))
2830 111 print *,"Environment variable ",string(:ilen(string))," not set."
2833 112 print *,"Error opening environment file!"
2835 call getenv(string,var)
2838 end subroutine mygetenv
2839 !-----------------------------------------------------------------------------
2841 !-----------------------------------------------------------------------------
2842 subroutine read_general_data(*)
2844 use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
2845 scelemode,TUBEmode,tor_mode
2847 use energy_data, only:distchainmax,tubeR0,tubecenter
2848 use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
2849 bordtubebot,tubebufthick,buftubebot,buftubetop
2851 ! include "DIMENSIONS"
2852 ! include "DIMENSIONS.ZSCOPT"
2853 ! include "DIMENSIONS.FREE"
2854 ! include "COMMON.TORSION"
2855 ! include "COMMON.INTERACT"
2856 ! include "COMMON.IOUNITS"
2857 ! include "COMMON.TIME1"
2858 ! include "COMMON.PROT"
2859 ! include "COMMON.PROTFILES"
2860 ! include "COMMON.CHAIN"
2861 ! include "COMMON.NAMES"
2862 ! include "COMMON.FFIELD"
2863 ! include "COMMON.ENEPS"
2864 ! include "COMMON.WEIGHTS"
2865 ! include "COMMON.FREE"
2866 ! include "COMMON.CONTROL"
2867 ! include "COMMON.ENERGIES"
2868 character(len=800) :: controlcard
2869 integer :: i,j,k,ii,n_ene_found
2870 integer :: ind,itype1,itype2,itypf,itypsc,itypp
2873 !el character*16 ucase
2874 character(len=16) :: key
2876 call card_concat(controlcard,.true.)
2877 call readi(controlcard,"N_ENE",n_eneW,max_eneW)
2878 if (n_eneW.gt.max_eneW) then
2879 write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
2883 call readi(controlcard,"NPARMSET",nparmset,1)
2884 !elwrite(iout,*)"in read_gen data"
2885 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
2886 call readi(controlcard,"IPARMPRINT",iparmprint,1)
2887 write (iout,*) "PARMPRINT",iparmprint
2888 if (nparmset.gt.max_parm) then
2889 write (iout,*) "Error: parameter out of range: NPARMSET",&
2893 !elwrite(iout,*)"in read_gen data"
2894 call readi(controlcard,"MAXIT",maxit,5000)
2895 call reada(controlcard,"FIMIN",fimin,1.0d-3)
2896 call readi(controlcard,"ENSEMBLES",ensembles,0)
2897 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
2898 write (iout,*) "Number of energy parameter sets",nparmset
2899 allocate(isampl(nparmset))
2900 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
2901 write (iout,*) "MaxSlice",MaxSlice
2902 call readi(controlcard,"NSLICE",nslice,1)
2903 !elwrite(iout,*)"in read_gen data"
2905 if (nslice.gt.MaxSlice) then
2906 write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
2910 write (iout,*) "Frequency of storing conformations",&
2911 (isampl(i),i=1,nparmset)
2912 write (iout,*) "Maxit",maxit," Fimin",fimin
2913 call readi(controlcard,"NQ",nQ,1)
2914 if (nQ.gt.MaxQ) then
2915 write (iout,*) "Error: parameter out of range: NQ",nq,&
2920 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
2921 call reada(controlcard,"DELTA",delta,1.0d-2)
2922 call readi(controlcard,"EINICHECK",einicheck,2)
2923 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
2924 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
2925 call readi(controlcard,"RESCALE",rescale_modeW,1)
2926 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2927 call reada(controlcard,'BOXY',boxysize,100.0d0)
2928 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2929 call readi(controlcard,"SCELEMODE",scelemode,0)
2930 print *,"SCELE",scelemode
2931 call readi(controlcard,'TORMODE',tor_mode,0)
2932 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2933 write(iout,*) "torsional and valence angle mode",tor_mode
2935 call readi(controlcard,'TUBEMOD',tubemode,0)
2937 if (TUBEmode.gt.0) then
2938 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
2939 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
2940 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
2941 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
2942 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
2943 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
2944 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
2945 buftubebot=bordtubebot+tubebufthick
2946 buftubetop=bordtubetop-tubebufthick
2948 ions=index(controlcard,"IONS").gt.0
2949 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2950 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2951 write(iout,*) "R_CUT_ELE=",r_cut_ele
2952 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
2953 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
2954 call readi(controlcard,'SYM',symetr,1)
2955 write (iout,*) "DISTCHAINMAX",distchainmax
2956 write (iout,*) "delta",delta
2957 write (iout,*) "einicheck",einicheck
2958 write (iout,*) "rescale_mode",rescale_modeW
2960 bxfile=index(controlcard,"BXFILE").gt.0
2961 cxfile=index(controlcard,"CXFILE").gt.0
2962 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
2964 histfile=index(controlcard,"HISTFILE").gt.0
2965 histout=index(controlcard,"HISTOUT").gt.0
2966 entfile=index(controlcard,"ENTFILE").gt.0
2967 zscfile=index(controlcard,"ZSCFILE").gt.0
2968 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
2969 write (iout,*) "with_dihed_constr ",with_dihed_constr
2970 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2972 end subroutine read_general_data
2973 !------------------------------------------------------------------------------
2974 subroutine read_efree(*)
2976 ! Read molecular data
2979 ! include 'DIMENSIONS'
2980 ! include 'DIMENSIONS.ZSCOPT'
2981 ! include 'DIMENSIONS.COMPAR'
2982 ! include 'DIMENSIONS.FREE'
2983 ! include 'COMMON.IOUNITS'
2984 ! include 'COMMON.TIME1'
2985 ! include 'COMMON.SBRIDGE'
2986 ! include 'COMMON.CONTROL'
2987 ! include 'COMMON.CHAIN'
2988 ! include 'COMMON.HEADER'
2989 ! include 'COMMON.GEO'
2990 ! include 'COMMON.FREE'
2991 character(len=320) :: controlcard !,ucase
2992 integer :: iparm,ib,i,j,npars
3002 ! call alloc_wham_arrays
3003 ! allocate(nT_h(nParmSet))
3004 ! allocate(replica(nParmSet))
3005 ! allocate(umbrella(nParmSet))
3006 ! allocate(read_iset(nParmSet))
3007 ! allocate(nT_h(nParmSet))
3011 call card_concat(controlcard,.true.)
3012 call readi(controlcard,'NT',nT_h(iparm),1)
3013 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
3015 if (nT_h(iparm).gt.MaxT_h) then
3016 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
3020 replica(iparm)=index(controlcard,"REPLICA").gt.0
3021 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
3022 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
3023 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
3024 replica(iparm)," umbrella ",umbrella(iparm),&
3025 " read_iset",read_iset(iparm)
3028 call card_concat(controlcard,.true.)
3029 call readi(controlcard,'NR',nR(ib,iparm),1)
3030 if (umbrella(iparm)) then
3033 nRR(ib,iparm)=nR(ib,iparm)
3035 if (nR(ib,iparm).gt.MaxR) then
3036 write (iout,*) "Error: parameter out of range: NR",&
3040 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
3041 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
3042 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
3045 call card_concat(controlcard,.true.)
3046 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
3048 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
3053 write (iout,*) "ib",ib," beta_h",&
3054 1.0d0/(0.001987*beta_h(ib,iparm))
3055 write (iout,*) "nR",nR(ib,iparm)
3056 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
3058 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
3059 "q0",(q0(j,i,ib,iparm),j=1,nQ)
3071 nR(ib,iparm)=nR(ib,1)
3072 if (umbrella(iparm)) then
3075 nRR(ib,iparm)=nR(ib,1)
3077 beta_h(ib,iparm)=beta_h(ib,1)
3079 f(i,ib,iparm)=f(i,ib,1)
3081 KH(j,i,ib,iparm)=KH(j,i,ib,1)
3082 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
3085 replica(iparm)=replica(1)
3086 umbrella(iparm)=umbrella(1)
3087 read_iset(iparm)=read_iset(1)
3094 end subroutine read_efree
3095 !-----------------------------------------------------------------------------
3096 subroutine read_protein_data(*)
3098 ! include "DIMENSIONS"
3099 ! include "DIMENSIONS.ZSCOPT"
3100 ! include "DIMENSIONS.FREE"
3104 integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
3105 ! include "COMMON.MPI"
3107 ! include "COMMON.CHAIN"
3108 ! include "COMMON.IOUNITS"
3109 ! include "COMMON.PROT"
3110 ! include "COMMON.PROTFILES"
3111 ! include "COMMON.NAMES"
3112 ! include "COMMON.FREE"
3113 ! include "COMMON.OBCINKA"
3114 character(len=64) :: nazwa
3115 character(len=16000) :: controlcard
3116 integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
3117 !el external ilen,iroof
3126 ! Read names of files with conformation data.
3127 if (replica(iparm)) then
3133 do ii=1,nRR(ib,iparm)
3134 write (iout,*) "Parameter set",iparm," temperature",ib,&
3137 call card_concat(controlcard,.true.)
3138 write (iout,*) controlcard(:ilen(controlcard))
3139 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
3140 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
3141 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
3142 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
3143 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
3144 maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
3145 call reada(controlcard,"TIME_START",&
3146 time_start_collect(ii,ib,iparm),0.0d0)
3147 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
3149 write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
3150 " rec_end",rec_end(ii,ib,iparm)
3151 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
3152 " time_end",time_end_collect(ii,ib,iparm)
3154 if (replica(iparm)) then
3155 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
3156 write (iout,*) "Number of trajectories",totraj(ii,iparm)
3159 if (nfile_bin(ii,ib,iparm).lt.2 &
3160 .and. nfile_asc(ii,ib,iparm).eq.0 &
3161 .and. nfile_cx(ii,ib,iparm).eq.0) then
3162 write (iout,*) "Error - no action specified!"
3165 if (nfile_bin(ii,ib,iparm).gt.0) then
3166 call card_concat(controlcard,.false.)
3167 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
3168 maxfile_prot,nfile_bin(ii,ib,iparm))
3170 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
3171 write(iout,*) (protfiles(i,1,ii,ib,iparm),&
3172 i=1,nfile_bin(ii,ib,iparm))
3175 if (nfile_asc(ii,ib,iparm).gt.0) then
3176 call card_concat(controlcard,.false.)
3177 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3178 maxfile_prot,nfile_asc(ii,ib,iparm))
3180 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
3181 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3182 i=1,nfile_asc(ii,ib,iparm))
3184 else if (nfile_cx(ii,ib,iparm).gt.0) then
3185 call card_concat(controlcard,.false.)
3186 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3187 maxfile_prot,nfile_cx(ii,ib,iparm))
3189 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
3190 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3191 i=1,nfile_cx(ii,ib,iparm))
3200 end subroutine read_protein_data
3201 !-------------------------------------------------------------------------------
3202 subroutine readsss(rekord,lancuch,wartosc,default)
3204 character*(*) :: rekord,lancuch,wartosc,default
3205 character(len=80) :: aux
3206 integer :: lenlan,lenrec,iread,ireade
3210 lenlan=ilen(lancuch)
3212 iread=index(rekord,lancuch(:lenlan)//"=")
3213 ! print *,"rekord",rekord," lancuch",lancuch
3214 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3215 if (iread.eq.0) then
3219 iread=iread+lenlan+1
3220 ! print *,"iread",iread
3221 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3222 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3224 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3226 ! print *,"iread",iread
3227 if (iread.gt.lenrec) then
3232 ! print *,"ireade",ireade
3233 do while (ireade.lt.lenrec .and. &
3234 .not.iblnk(rekord(ireade:ireade)))
3237 wartosc=rekord(iread:ireade)
3239 end subroutine readsss
3240 !----------------------------------------------------------------------------
3241 subroutine multreads(rekord,lancuch,tablica,dim,default)
3244 character*(*) rekord,lancuch,tablica(dim),default
3245 character(len=80) :: aux
3246 integer :: lenlan,lenrec,iread,ireade
3253 lenlan=ilen(lancuch)
3255 iread=index(rekord,lancuch(:lenlan)//"=")
3256 ! print *,"rekord",rekord," lancuch",lancuch
3257 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3258 if (iread.eq.0) return
3259 iread=iread+lenlan+1
3261 ! print *,"iread",iread
3262 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3263 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3265 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3267 ! print *,"iread",iread
3268 if (iread.gt.lenrec) return
3270 ! print *,"ireade",ireade
3271 do while (ireade.lt.lenrec .and. &
3272 .not.iblnk(rekord(ireade:ireade)))
3275 tablica(i)=rekord(iread:ireade)
3278 end subroutine multreads
3279 !----------------------------------------------------------------------------
3280 subroutine split_string(rekord,tablica,dim,nsub)
3282 integer :: dim,nsub,i,ii,ll,kk
3283 character*(*) tablica(dim)
3284 character*(*) rekord
3294 ! Find the start of term name
3296 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
3299 ! Parse the name into TABLICA(i) until blank found
3300 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
3302 tablica(i)(kk:kk)=rekord(ii:ii)
3305 if (kk.gt.0) nsub=nsub+1
3306 if (ii.gt.ll) return
3309 end subroutine split_string
3310 !--------------------------------------------------------------------------------
3312 !--------------------------------------------------------------------------------
3313 subroutine read_compar
3315 ! Read molecular data
3317 use conform_compar, only:alloc_compar_arrays
3318 use control_data, only:pdbref
3319 use geometry_data, only:deg2rad,rad2deg
3321 ! include 'DIMENSIONS'
3322 ! include 'DIMENSIONS.ZSCOPT'
3323 ! include 'DIMENSIONS.COMPAR'
3324 ! include 'DIMENSIONS.FREE'
3325 ! include 'COMMON.IOUNITS'
3326 ! include 'COMMON.TIME1'
3327 ! include 'COMMON.SBRIDGE'
3328 ! include 'COMMON.CONTROL'
3329 ! include 'COMMON.COMPAR'
3330 ! include 'COMMON.CHAIN'
3331 ! include 'COMMON.HEADER'
3332 ! include 'COMMON.GEO'
3333 ! include 'COMMON.FREE'
3334 character(len=320) :: controlcard !,ucase
3335 character(len=64) :: wfile
3339 !elwrite(iout,*)"jestesmy w read_compar"
3340 call card_concat(controlcard,.true.)
3341 pdbref=(index(controlcard,'PDBREF').gt.0)
3342 call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
3343 call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
3344 call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
3345 call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
3346 verbose = index(controlcard,"VERBOSE").gt.0
3347 lgrp=index(controlcard,"STATIN").gt.0
3348 lgrp_out=index(controlcard,"STATOUT").gt.0
3349 merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
3350 binary = index(controlcard,"BINARY").gt.0
3351 rmscut_base_up=rmscut_base_up/50
3352 rmscut_base_low=rmscut_base_low/50
3353 call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
3354 call readi(controlcard,'NLEVEL',nlevel,1)
3355 if (nlevel.lt.0) then
3357 call alloc_compar_arrays(maxfrag,1)
3360 allocate(nfrag(nlevel))
3362 ! Read the data pertaining to elementary fragments (level 1)
3363 call readi(controlcard,'NFRAG',nfrag(1),0)
3364 write(iout,*)"nfrag(1)",nfrag(1)
3365 call alloc_compar_arrays(nfrag(1),nlevel)
3367 call card_concat(controlcard,.true.)
3368 write (iout,*) controlcard(:ilen(controlcard))
3369 call readi(controlcard,'NPIECE',npiece(j,1),0)
3370 call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
3371 call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
3372 call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
3373 call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
3374 call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
3375 call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
3376 call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
3377 call readi(controlcard,'RMS',irms(j,1),0)
3378 call readi(controlcard,'LOCAL',iloc(j),1)
3379 call readi(controlcard,'ELCONT',ielecont(j,1),1)
3380 if (ielecont(j,1).eq.0) then
3381 call readi(controlcard,'SCCONT',isccont(j,1),1)
3383 ang_cut(j)=ang_cut(j)*deg2rad
3384 ang_cut1(j)=ang_cut1(j)*deg2rad
3386 call card_concat(controlcard,.true.)
3387 call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
3388 call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
3390 write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
3391 (ifrag(1,k,j),ifrag(2,k,j),&
3392 k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
3393 " ang_cut1",ang_cut1(j)*rad2deg
3394 write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
3395 write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
3396 write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
3397 " ilocal",iloc(j)," isccont",isccont(j,1)
3399 ! Read data pertaning to higher levels
3401 call card_concat(controlcard,.true.)
3402 call readi(controlcard,'NFRAG',NFRAG(i),0)
3403 write (iout,*) "i",i," nfrag",nfrag(i)
3405 call card_concat(controlcard,.true.)
3407 call readi(controlcard,'ELCONT',ielecont(j,i),0)
3408 if (ielecont(j,i).eq.0) then
3409 call readi(controlcard,'SCCONT',isccont(j,i),1)
3411 call readi(controlcard,'RMS',irms(j,i),0)
3417 call readi(controlcard,'NPIECE',npiece(j,i),0)
3418 call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
3419 call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
3420 call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
3422 call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
3423 call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
3424 write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
3425 n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
3426 " isccont",isccont(j,i)," irms",irms(j,i)
3427 write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
3428 write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
3429 write(iout,*)"nc_frac",nc_fragm(j,i),&
3430 " nc_req",nc_req_setf(j,i)
3433 if (binary) write (iout,*) "Classes written in binary format."
3436 call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
3437 call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
3438 call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
3439 call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
3440 call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
3441 call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
3442 call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
3443 call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
3444 call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
3445 call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
3446 call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
3447 call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
3448 call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
3449 call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
3450 call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
3451 call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
3452 call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
3453 call readi(controlcard,'RMS_SINGLE',irms_single,0)
3454 call readi(controlcard,'CONT_SINGLE',icont_single,1)
3455 call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
3456 call readi(controlcard,'RMS_PAIR',irms_pair,0)
3457 call readi(controlcard,'CONT_PAIR',icont_pair,1)
3458 call readi(controlcard,'SPLIT_BET',isplit_bet,0)
3459 angcut_hel=angcut_hel*deg2rad
3460 angcut1_hel=angcut1_hel*deg2rad
3461 angcut_bet=angcut_bet*deg2rad
3462 angcut1_bet=angcut1_bet*deg2rad
3463 angcut_strand=angcut_strand*deg2rad
3464 angcut1_strand=angcut1_strand*deg2rad
3465 write (iout,*) "Automatic detection of structural elements"
3466 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
3467 ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
3468 ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
3469 ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
3470 ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
3471 ' SPLIT_BET',isplit_bet
3472 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
3473 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
3474 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
3475 ' MAXANG_HEL',angcut1_hel*rad2deg
3476 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
3477 ' MAXANG_BET',angcut1_bet*rad2deg
3478 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
3479 ' MAXANG_STRAND',angcut1_strand*rad2deg
3480 write (iout,*) 'FRAC_MIN',frac_min_set
3482 end subroutine read_compar
3483 !--------------------------------------------------------------------------------
3485 !--------------------------------------------------------------------------------
3486 subroutine read_ref_structure(*)
3488 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
3491 use control_data, only:pdbref
3492 use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
3493 nstart_sup,nstart_seq,nperm,nres0
3494 use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
3495 use compare, only:seq_comp !,contact,elecont
3496 use geometry, only:chainbuild,dist
3497 use io_config, only:readpdb
3499 use conform_compar, only:contact,elecont
3501 ! include 'DIMENSIONS'
3502 ! include 'DIMENSIONS.ZSCOPT'
3503 ! include 'DIMENSIONS.COMPAR'
3504 ! include 'COMMON.IOUNITS'
3505 ! include 'COMMON.GEO'
3506 ! include 'COMMON.VAR'
3507 ! include 'COMMON.INTERACT'
3508 ! include 'COMMON.LOCAL'
3509 ! include 'COMMON.NAMES'
3510 ! include 'COMMON.CHAIN'
3511 ! include 'COMMON.FFIELD'
3512 ! include 'COMMON.SBRIDGE'
3513 ! include 'COMMON.HEADER'
3514 ! include 'COMMON.CONTROL'
3515 ! include 'COMMON.CONTACTS1'
3516 ! include 'COMMON.PEPTCONT'
3517 ! include 'COMMON.TIME1'
3518 ! include 'COMMON.COMPAR'
3519 character(len=4) :: sequence(nres)
3521 !el real(kind=8) :: x(maxvar)
3522 integer :: itype_pdb(nres,5)
3523 !el logical seq_comp
3524 integer :: i,j,k,nres_pdb,iaux,mnum
3525 real(kind=8) :: ddsc !el,dist
3526 integer :: kkk !,ilen
3530 write (iout,*) "pdbref",pdbref
3532 read(inp,'(a)') pdbfile
3533 write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
3534 pdbfile(:ilen(pdbfile))
3535 open(ipdbin,file=pdbfile,status='old',err=33)
3537 33 write (iout,'(a)') 'Error opening PDB file.'
3542 itype_pdb(i,mnum)=itype(i,mnum)
3548 iaux=itype_pdb(i,mnum)
3549 itype_pdb(i,mnum)=itype(i,mnum)
3557 if (nsup.le.(nct-nnt+1)) then
3558 do i=0,nct-nnt+1-nsup
3559 if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
3561 do j=nnt+nsup-1,nnt,-1
3563 cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
3566 do j=nnt+nsup-1,nnt,-1
3568 cref(k,j+i,kkk)=cref(k,j,kkk)
3570 write(iout,*) "J",j,"J+I",j+i
3571 phi_ref(j+i)=phi_ref(j)
3572 theta_ref(j+i)=theta_ref(j)
3573 alph_ref(j+i)=alph_ref(j)
3574 omeg_ref(j+i)=omeg_ref(j)
3578 write (iout,'(i5,3f10.5,5x,3f10.5)') &
3579 j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
3587 write (iout,'(a)') &
3588 'Error - sequences to be superposed do not match.'
3591 do i=0,nsup-(nct-nnt+1)
3592 if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
3595 nstart_sup=nstart_sup+i
3600 write (iout,'(a)') &
3601 'Error - sequences to be superposed do not match.'
3605 write (iout,'(a,i5)') &
3606 'Experimental structure begins at residue',nstart_seq
3608 call read_angles(inp,*38)
3610 38 write (iout,'(a)') 'Error reading reference structure.'
3619 cref(j,i,kkk)=c(j,i)
3623 nend_sup=nstart_sup+nsup-1
3626 c(j,i)=cref(j,i,kkk)
3632 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
3634 if (itype(i,mnum).ne.10) then
3635 ddsc = dist(i,nres+i)
3637 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
3641 dc_norm(j,nres+i)=0.0d0
3644 ! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
3645 ! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
3646 ! dc_norm(3,nres+i)**2
3648 dc(j,i)=c(j,i+1)-c(j,i)
3652 dc_norm(j,i)=dc(j,i)/ddsc
3655 ! print *,"Calling contact"
3656 call contact(.true.,ncont_ref,icont_ref(1,1),&
3657 nstart_sup,nend_sup)
3658 ! print *,"Calling elecont"
3659 call elecont(.true.,ncont_pept_ref,&
3660 icont_pept_ref(1,1),&
3661 nstart_sup,nend_sup)
3662 write (iout,'(a,i3,a,i3,a,i3,a)') &
3663 'Number of residues to be superposed:',nsup,&
3664 ' (from residue',nstart_sup,' to residue',&
3667 end subroutine read_ref_structure
3668 !--------------------------------------------------------------------------------
3670 !--------------------------------------------------------------------------------
3671 subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
3673 use geometry_data, only:nres,c
3674 use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
3675 ! implicit real*8 (a-h,o-z)
3676 ! include 'DIMENSIONS'
3677 ! include 'DIMENSIONS.ZSCOPT'
3678 ! include 'COMMON.CHAIN'
3679 ! include 'COMMON.INTERACT'
3680 ! include 'COMMON.NAMES'
3681 ! include 'COMMON.IOUNITS'
3682 ! include 'COMMON.HEADER'
3683 ! include 'COMMON.SBRIDGE'
3684 character(len=50) :: tytul
3685 character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
3686 'D','E','F','G','H','I','J'/),shape(chainid))
3687 integer,dimension(nres) :: ica !(maxres)
3688 real(kind=8) :: temp,efree,etot,entropy,rmsdev
3689 integer :: ii,i,j,iti,ires,iatom,ichain,mnum
3690 write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
3692 write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
3694 write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
3702 if (iti.eq.ntyp1) then
3705 write (ipdb,'(a)') 'TER'
3710 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
3714 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
3715 ires,(c(j,nres+i),j=1,3)
3719 write (ipdb,'(a)') 'TER'
3722 if (itype(i,mnum).eq.ntyp1) cycle
3723 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3724 write (ipdb,30) ica(i),ica(i+1)
3725 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3726 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
3727 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
3728 write (ipdb,30) ica(i),ica(i)+1
3731 if (itype(nct,molnum(nct)).ne.10) then
3732 write (ipdb,30) ica(nct),ica(nct)+1
3735 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
3737 write (ipdb,'(a6)') 'ENDMDL'
3738 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3739 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3740 30 FORMAT ('CONECT',8I5)
3742 end subroutine pdboutW
3744 !------------------------------------------------------------------------------
3746 !-----------------------------------------------------------------------------
3747 !-----------------------------------------------------------------------------