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
487 ! print *,"KURWA",ww(48)
488 ! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO "
489 ! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",&
490 ! "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",&
491 ! "WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",&
492 ! "WCATPROT ","WCATCAT
493 ! +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
494 ! +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
495 ! +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
496 ! +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
501 allocate(weights(n_ene))
516 weights(15)=wstrain !0
519 weights(18)=0 !scal14 !
523 weights(26)= wvdwpp_nucl
530 weights(32) =wbond_nucl
531 weights(33) =wang_nucl
533 weights(35) =wtor_nucl
534 weights(36) =wtor_d_nucl
535 weights(37) =wcorr_nucl
536 weights(38) =wcorr3_nucl
538 weights(42) =wcatprot
543 call card_concat(controlcard,.false.)
545 ! Return if not own parameters
547 if (iparm.ne.myparm .and. separate_parset) return
549 call reads(controlcard,"BONDPAR",bondname_t,bondname)
550 open (ibond,file=bondname_t,status='old')
552 call reads(controlcard,"THETPAR",thetname_t,thetname)
553 open (ithep,file=thetname_t,status='old')
555 call reads(controlcard,"ROTPAR",rotname_t,rotname)
556 open (irotam,file=rotname_t,status='old')
558 call reads(controlcard,"TORPAR",torname_t,torname)
559 open (itorp,file=torname_t,status='old')
561 call reads(controlcard,"TORDPAR",tordname_t,tordname)
562 open (itordp,file=tordname_t,status='old')
564 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
565 open (isccor,file=sccorname_t,status='old')
567 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
568 open (ifourier,file=fouriername_t,status='old')
570 call reads(controlcard,"ELEPAR",elename_t,elename)
571 open (ielep,file=elename_t,status='old')
573 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
574 open (isidep,file=sidename_t,status='old')
576 call reads(controlcard,"SCPPAR",scpname_t,scpname)
577 open (iscpp,file=scpname_t,status='old')
579 call getenv_loc('IONPAR',ionname)
580 open (iion,file=ionname,status='old')
582 write (iout,*) "Parameter set:",iparm
583 write (iout,*) "Energy-term weights:"
585 write (iout,'(a16,f10.5)') wname(i),ww(i)
587 write (iout,*) "Sidechain potential file : ",&
588 sidename_t(:ilen(sidename_t))
590 write (iout,*) "SCp potential file : ",&
591 scpname_t(:ilen(scpname_t))
593 write (iout,*) "Electrostatic potential file : ",&
594 elename_t(:ilen(elename_t))
595 write (iout,*) "Cumulant coefficient file : ",&
596 fouriername_t(:ilen(fouriername_t))
597 write (iout,*) "Torsional parameter file : ",&
598 torname_t(:ilen(torname_t))
599 write (iout,*) "Double torsional parameter file : ",&
600 tordname_t(:ilen(tordname_t))
601 write (iout,*) "Backbone-rotamer parameter file : ",&
602 sccorname(:ilen(sccorname))
603 write (iout,*) "Bond & inertia constant file : ",&
604 bondname_t(:ilen(bondname_t))
605 write (iout,*) "Bending parameter file : ",&
606 thetname_t(:ilen(thetname_t))
607 write (iout,*) "Rotamer parameter file : ",&
608 rotname_t(:ilen(rotname_t))
611 ! Read the virtual-bond parameters, masses, and moments of inertia
612 ! and Stokes' radii of the peptide group and side chains
614 allocate(dsc(ntyp1)) !(ntyp1)
615 allocate(dsc_inv(ntyp1)) !(ntyp1)
616 allocate(nbondterm(ntyp)) !(ntyp)
617 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
618 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
619 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
620 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
621 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
623 !el allocate(msc(ntyp+1)) !(ntyp+1)
624 !el allocate(isc(ntyp+1)) !(ntyp+1)
625 !el allocate(restok(ntyp+1)) !(ntyp+1)
626 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
629 read (ibond,*) vbldp0,akp
632 read (ibond,*) vbldsc0(1,i),aksc(1,i)
633 dsc(i) = vbldsc0(1,i)
637 dsc_inv(i)=1.0D0/dsc(i)
641 read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
643 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
645 dsc(i) = vbldsc0(1,i)
649 dsc_inv(i)=1.0D0/dsc(i)
654 write(iout,'(/a/)')"Force constants virtual bonds:"
655 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
657 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
659 write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
660 vbldsc0(1,i),aksc(1,i),abond0(1,i)
662 write (iout,'(13x,3f10.5)') &
663 vbldsc0(j,i),aksc(j,i),abond0(j,i)
667 if (.not. allocated(msc)) allocate(msc(ntyp1,5))
668 if (.not. allocated(restok)) allocate(restok(ntyp1,5))
671 read(iion,*) msc(i,5),restok(i,5)
672 print *,msc(i,5),restok(i,5)
676 read (iion,*) ncatprotparm
677 allocate(catprm(ncatprotparm,4))
679 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
682 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
685 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2)
686 ! dsc(i) = vbldsc0_nucl(1,i)
690 ! dsc_inv(i)=1.0D0/dsc(i)
695 !----------------------------------------------------
696 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
697 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
698 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
699 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
700 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
701 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
707 athet(j,i,ichir1,ichir2)=0.0D0
708 bthet(j,i,ichir1,ichir2)=0.0D0
722 !elwrite(iout,*) "parmread kontrol"
726 ! Read the parameters of the probability distribution/energy expression
727 ! of the virtual-bond valence angles theta
730 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
731 (bthet(j,i,1,1),j=1,2)
732 read (ithep,*) (polthet(j,i),j=0,3)
733 !elwrite(iout,*) "parmread kontrol in cryst_theta"
734 read (ithep,*) (gthet(j,i),j=1,3)
735 !elwrite(iout,*) "parmread kontrol in cryst_theta"
736 read (ithep,*) theta0(i),sig0(i),sigc0(i)
738 !elwrite(iout,*) "parmread kontrol in cryst_theta"
740 !elwrite(iout,*) "parmread kontrol in cryst_theta"
742 athet(1,i,1,-1)=athet(1,i,1,1)
743 athet(2,i,1,-1)=athet(2,i,1,1)
744 bthet(1,i,1,-1)=-bthet(1,i,1,1)
745 bthet(2,i,1,-1)=-bthet(2,i,1,1)
746 athet(1,i,-1,1)=-athet(1,i,1,1)
747 athet(2,i,-1,1)=-athet(2,i,1,1)
748 bthet(1,i,-1,1)=bthet(1,i,1,1)
749 bthet(2,i,-1,1)=bthet(2,i,1,1)
751 !elwrite(iout,*) "parmread kontrol in cryst_theta"
754 athet(1,i,-1,-1)=athet(1,-i,1,1)
755 athet(2,i,-1,-1)=-athet(2,-i,1,1)
756 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
757 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
758 athet(1,i,-1,1)=athet(1,-i,1,1)
759 athet(2,i,-1,1)=-athet(2,-i,1,1)
760 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
761 bthet(2,i,-1,1)=bthet(2,-i,1,1)
762 athet(1,i,1,-1)=-athet(1,-i,1,1)
763 athet(2,i,1,-1)=athet(2,-i,1,1)
764 bthet(1,i,1,-1)=bthet(1,-i,1,1)
765 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
770 polthet(j,i)=polthet(j,-i)
773 gthet(j,i)=gthet(j,-i)
776 !elwrite(iout,*) "parmread kontrol in cryst_theta"
778 !elwrite(iout,*) "parmread kontrol in cryst_theta"
781 ! & 'Parameters of the virtual-bond valence angles:'
782 ! write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
783 ! & ' ATHETA0 ',' A1 ',' A2 ',
786 ! write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
787 ! & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
789 ! write (iout,'(/a/9x,5a/79(1h-))')
790 ! & 'Parameters of the expression for sigma(theta_c):',
791 ! & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
792 ! & ' ALPH3 ',' SIGMA0C '
794 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
795 ! & (polthet(j,i),j=0,3),sigc0(i)
797 ! write (iout,'(/a/9x,5a/79(1h-))')
798 ! & 'Parameters of the second gaussian:',
799 ! & ' THETA0 ',' SIGMA0 ',' G1 ',
802 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
803 ! & sig0(i),(gthet(j,i),j=1,3)
806 'Parameters of the virtual-bond valence angles:'
807 write (iout,'(/a/9x,5a/79(1h-))') &
808 'Coefficients of expansion',&
809 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
810 ' b1*10^1 ',' b2*10^1 '
812 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
813 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
814 (10*bthet(j,i,1,1),j=1,2)
816 write (iout,'(/a/9x,5a/79(1h-))') &
817 'Parameters of the expression for sigma(theta_c):',&
818 ' alpha0 ',' alph1 ',' alph2 ',&
819 ' alhp3 ',' sigma0c '
821 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
822 (polthet(j,i),j=0,3),sigc0(i)
824 write (iout,'(/a/9x,5a/79(1h-))') &
825 'Parameters of the second gaussian:',&
826 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
829 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
830 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
835 ! Read the parameters of Utheta determined from ab initio surfaces
836 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
838 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
839 ! write (iout,*) "tu dochodze"
840 IF (tor_mode.eq.0) THEN
841 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
842 read (ithep,*) nthetyp,ntheterm,ntheterm2,&
843 ntheterm3,nsingle,ndouble
844 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
846 !----------------------------------------------------
847 ! allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
848 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
849 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
850 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
851 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
852 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
853 !(maxtheterm,-maxthetyp1:maxthetyp1,&
854 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
855 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
856 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
857 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
858 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
859 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
860 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
861 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
862 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
863 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
864 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
865 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
866 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
867 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
868 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
869 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
870 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
873 read (ithep,*) (ithetyp(i),i=1,ntyp1)
875 ithetyp(i)=-ithetyp(-i)
877 ! write (iout,*) "tu dochodze"
878 aa0thet(:,:,:,:)=0.0d0
879 aathet(:,:,:,:,:)=0.0d0
880 bbthet(:,:,:,:,:,:)=0.0d0
881 ccthet(:,:,:,:,:,:)=0.0d0
882 ddthet(:,:,:,:,:,:)=0.0d0
883 eethet(:,:,:,:,:,:)=0.0d0
884 ffthet(:,:,:,:,:,:,:)=0.0d0
885 ggthet(:,:,:,:,:,:,:)=0.0d0
889 do j=-nthetyp,nthetyp
890 do k=-nthetyp,nthetyp
891 read (ithep,'(6a)') res1
892 read (ithep,*) aa0thet(i,j,k,iblock)
893 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
895 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
896 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
897 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
898 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
901 (((ffthet(llll,lll,ll,i,j,k,iblock),&
902 ffthet(lll,llll,ll,i,j,k,iblock),&
903 ggthet(llll,lll,ll,i,j,k,iblock),&
904 ggthet(lll,llll,ll,i,j,k,iblock),&
905 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
910 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
911 ! coefficients of theta-and-gamma-dependent terms are zero.
916 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
917 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
919 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
920 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
923 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
925 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
928 ! Substitution for D aminoacids from symmetry.
931 do j=-nthetyp,nthetyp
932 do k=-nthetyp,nthetyp
933 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
935 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
939 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
940 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
941 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
942 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
948 ffthet(llll,lll,ll,i,j,k,iblock)= &
949 ffthet(llll,lll,ll,-i,-j,-k,iblock)
950 ffthet(lll,llll,ll,i,j,k,iblock)= &
951 ffthet(lll,llll,ll,-i,-j,-k,iblock)
952 ggthet(llll,lll,ll,i,j,k,iblock)= &
953 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
954 ggthet(lll,llll,ll,i,j,k,iblock)= &
955 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
965 ! Control printout of the coefficients of virtual-bond-angle potentials
969 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
973 write (iout,'(//4a)') &
974 'Type ',onelett(i),onelett(j),onelett(k)
975 write (iout,'(//a,10x,a)') " l","a[l]"
976 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
977 write (iout,'(i2,1pe15.5)') &
978 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
980 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
981 "b",l,"c",l,"d",l,"e",l
983 write (iout,'(i2,4(1pe15.5))') m,&
984 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
985 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
989 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
990 "f+",l,"f-",l,"g+",l,"g-",l
993 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
994 ffthet(n,m,l,i,j,k,iblock),&
995 ffthet(m,n,l,i,j,k,iblock),&
996 ggthet(n,m,l,i,j,k,iblock),&
997 ggthet(m,n,l,i,j,k,iblock)
1008 !C here will be the apropriate recalibrating for D-aminoacid
1009 read (ithep,*) nthetyp
1010 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1011 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1012 do i=-nthetyp+1,nthetyp-1
1013 read (ithep,*) nbend_kcc_Tb(i)
1014 do j=0,nbend_kcc_Tb(i)
1015 read (ithep,*) ijunk,v1bend_chyb(j,i)
1019 write (iout,'(a)') &
1020 "Parameters of the valence-only potentials"
1021 do i=-nthetyp+1,nthetyp-1
1022 write (iout,'(2a)') "Type ",toronelet(i)
1023 do k=0,nbend_kcc_Tb(i)
1024 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1031 !--------------- Reading theta parameters for nucleic acid-------
1032 read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,&
1033 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1034 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1035 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1036 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1037 nthetyp_nucl+1,nthetyp_nucl+1))
1038 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1039 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1040 nthetyp_nucl+1,nthetyp_nucl+1))
1041 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1042 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1043 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1044 nthetyp_nucl+1,nthetyp_nucl+1))
1045 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1046 nthetyp_nucl+1,nthetyp_nucl+1))
1047 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1048 nthetyp_nucl+1,nthetyp_nucl+1))
1049 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1050 nthetyp_nucl+1,nthetyp_nucl+1))
1051 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1052 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1053 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1054 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1055 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1056 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1058 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1059 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1061 read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1063 aa0thet_nucl(:,:,:)=0.0d0
1064 aathet_nucl(:,:,:,:)=0.0d0
1065 bbthet_nucl(:,:,:,:,:)=0.0d0
1066 ccthet_nucl(:,:,:,:,:)=0.0d0
1067 ddthet_nucl(:,:,:,:,:)=0.0d0
1068 eethet_nucl(:,:,:,:,:)=0.0d0
1069 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1070 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1075 read (ithep_nucl,'(3a)') t1,t2,t3
1076 read (ithep_nucl,*) aa0thet_nucl(i,j,k)
1077 read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1078 read (ithep_nucl,*) &
1079 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1080 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1081 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1082 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1083 read (ithep_nucl,*) &
1084 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1085 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1086 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1092 !-------------------------------------------
1093 allocate(nlob(ntyp1)) !(ntyp1)
1094 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1095 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1096 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1114 gaussc(l,k,j,i)=0.0D0
1122 ! Read the parameters of the probability distribution/energy expression
1123 ! of the side chains.
1126 !c write (iout,*) "tu dochodze",i
1127 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
1131 dsc_inv(i)=1.0D0/dsc(i)
1142 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
1143 censc(1,1,-i)=censc(1,1,i)
1144 censc(2,1,-i)=censc(2,1,i)
1145 censc(3,1,-i)=-censc(3,1,i)
1147 read (irotam,*) bsc(j,i)
1148 read (irotam,*) (censc(k,j,i),k=1,3),&
1149 ((blower(k,l,j),l=1,k),k=1,3)
1150 censc(1,j,-i)=censc(1,j,i)
1151 censc(2,j,-i)=censc(2,j,i)
1152 censc(3,j,-i)=-censc(3,j,i)
1153 ! BSC is amplitude of Gaussian
1160 akl=akl+blower(k,m,j)*blower(l,m,j)
1164 if (((k.eq.3).and.(l.ne.3)) &
1165 .or.((l.eq.3).and.(k.ne.3))) then
1166 gaussc(k,l,j,-i)=-akl
1167 gaussc(l,k,j,-i)=-akl
1169 gaussc(k,l,j,-i)=akl
1170 gaussc(l,k,j,-i)=akl
1179 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1182 if (nlobi.gt.0) then
1183 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1184 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1185 ! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1186 ! write (iout,'(a,f10.4,4(16x,f10.4))')
1187 ! & 'Center ',(bsc(j,i),j=1,nlobi)
1188 ! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
1189 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1190 'log h',(bsc(j,i),j=1,nlobi)
1191 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1192 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1193 ! write (iout,'(a)')
1199 ! blower(k,l,j)=gaussc(ind,j,i)
1204 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1205 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1212 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1213 ! added by Urszula Kozlowska 07/11/2007
1215 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1223 read(irotam,*) sc_parmin(j,i)
1229 !---------reading nucleic acid parameters for rotamers-------------------
1230 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1231 do i=1,ntyp_molec(2)
1232 read (irotam_nucl,*)
1234 read(irotam_nucl,*) sc_parmin_nucl(j,i)
1240 write (iout,*) "Base rotamer parameters"
1241 do i=1,ntyp_molec(2)
1242 write (iout,'(a)') restyp(i,2)
1243 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1248 read (ifourier,*) nloctyp
1249 !el write(iout,*)"nloctyp",nloctyp
1250 SPLIT_FOURIERTOR = nloctyp.lt.0
1251 nloctyp = iabs(nloctyp)
1253 if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1))
1254 print *,"shape",shape(itype2loc)
1255 allocate(iloctyp(-nloctyp:nloctyp))
1256 allocate(bnew1(3,2,-nloctyp:nloctyp))
1257 allocate(bnew2(3,2,-nloctyp:nloctyp))
1258 allocate(ccnew(3,2,-nloctyp:nloctyp))
1259 allocate(ddnew(3,2,-nloctyp:nloctyp))
1260 allocate(e0new(3,-nloctyp:nloctyp))
1261 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1262 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1263 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1264 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1265 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1266 allocate(e0newtor(3,-nloctyp:nloctyp))
1267 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1269 read (ifourier,*) (itype2loc(i),i=1,ntyp)
1270 read (ifourier,*) (iloctyp(i),i=0,nloctyp-1)
1271 itype2loc(ntyp1)=nloctyp
1272 iloctyp(nloctyp)=ntyp1
1274 itype2loc(-i)=-itype2loc(i)
1277 allocate(iloctyp(-nloctyp:nloctyp))
1278 allocate(itype2loc(-ntyp1:ntyp1))
1285 iloctyp(-i)=-iloctyp(i)
1287 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1288 !c write (iout,*) "nloctyp",nloctyp,
1289 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1290 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1291 !c write (iout,*) "nloctyp",nloctyp,
1292 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1295 !c write (iout,*) "NEWCORR",i
1299 read (ifourier,*) bnew1(ii,j,i)
1302 !c write (iout,*) "NEWCORR BNEW1"
1303 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1306 read (ifourier,*) bnew2(ii,j,i)
1309 !c write (iout,*) "NEWCORR BNEW2"
1310 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1312 read (ifourier,*) ccnew(kk,1,i)
1313 read (ifourier,*) ccnew(kk,2,i)
1315 !c write (iout,*) "NEWCORR CCNEW"
1316 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1318 read (ifourier,*) ddnew(kk,1,i)
1319 read (ifourier,*) ddnew(kk,2,i)
1321 !c write (iout,*) "NEWCORR DDNEW"
1322 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1326 read (ifourier,*) eenew(ii,jj,kk,i)
1330 !c write (iout,*) "NEWCORR EENEW1"
1331 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1333 read (ifourier,*) e0new(ii,i)
1335 !c write (iout,*) (e0new(ii,i),ii=1,3)
1337 !c write (iout,*) "NEWCORR EENEW"
1340 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1341 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1342 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1343 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1348 bnew1(ii,1,-i)= bnew1(ii,1,i)
1349 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1350 bnew2(ii,1,-i)= bnew2(ii,1,i)
1351 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1354 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1355 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1356 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1357 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1358 ccnew(ii,1,-i)=ccnew(ii,1,i)
1359 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1360 ddnew(ii,1,-i)=ddnew(ii,1,i)
1361 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1363 e0new(1,-i)= e0new(1,i)
1364 e0new(2,-i)=-e0new(2,i)
1365 e0new(3,-i)=-e0new(3,i)
1367 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1368 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1369 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1370 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1374 write (iout,'(a)') "Coefficients of the multibody terms"
1375 do i=-nloctyp+1,nloctyp-1
1376 write (iout,*) "Type: ",onelet(iloctyp(i))
1377 write (iout,*) "Coefficients of the expansion of B1"
1379 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1381 write (iout,*) "Coefficients of the expansion of B2"
1383 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1385 write (iout,*) "Coefficients of the expansion of C"
1386 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1387 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1388 write (iout,*) "Coefficients of the expansion of D"
1389 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1390 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1391 write (iout,*) "Coefficients of the expansion of E"
1392 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1395 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1400 IF (SPLIT_FOURIERTOR) THEN
1402 !c write (iout,*) "NEWCORR TOR",i
1406 read (ifourier,*) bnew1tor(ii,j,i)
1409 !c write (iout,*) "NEWCORR BNEW1 TOR"
1410 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1413 read (ifourier,*) bnew2tor(ii,j,i)
1416 !c write (iout,*) "NEWCORR BNEW2 TOR"
1417 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1419 read (ifourier,*) ccnewtor(kk,1,i)
1420 read (ifourier,*) ccnewtor(kk,2,i)
1422 !c write (iout,*) "NEWCORR CCNEW TOR"
1423 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1425 read (ifourier,*) ddnewtor(kk,1,i)
1426 read (ifourier,*) ddnewtor(kk,2,i)
1428 !c write (iout,*) "NEWCORR DDNEW TOR"
1429 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1433 read (ifourier,*) eenewtor(ii,jj,kk,i)
1437 !c write (iout,*) "NEWCORR EENEW1 TOR"
1438 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1440 read (ifourier,*) e0newtor(ii,i)
1442 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1444 !c write (iout,*) "NEWCORR EENEW TOR"
1447 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1448 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1449 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1450 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1455 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1456 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1457 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1458 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1461 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1462 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1463 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1464 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1466 e0newtor(1,-i)= e0newtor(1,i)
1467 e0newtor(2,-i)=-e0newtor(2,i)
1468 e0newtor(3,-i)=-e0newtor(3,i)
1470 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1471 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1472 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1473 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1477 write (iout,'(a)') &
1478 "Single-body coefficients of the torsional potentials"
1479 do i=-nloctyp+1,nloctyp-1
1480 write (iout,*) "Type: ",onelet(iloctyp(i))
1481 write (iout,*) "Coefficients of the expansion of B1tor"
1483 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1484 j,(bnew1tor(k,j,i),k=1,3)
1486 write (iout,*) "Coefficients of the expansion of B2tor"
1488 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1489 j,(bnew2tor(k,j,i),k=1,3)
1491 write (iout,*) "Coefficients of the expansion of Ctor"
1492 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1493 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1494 write (iout,*) "Coefficients of the expansion of Dtor"
1495 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1496 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1497 write (iout,*) "Coefficients of the expansion of Etor"
1498 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1501 write (iout,'(1hE,2i1,2f10.5)') &
1502 j,k,(eenewtor(l,j,k,i),l=1,2)
1508 do i=-nloctyp+1,nloctyp-1
1511 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1516 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1520 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1521 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1522 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1523 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1526 ENDIF !SPLIT_FOURIER_TOR
1528 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1529 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1530 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1531 allocate(b(13,-nloctyp-1:nloctyp+1))
1533 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1536 read (ifourier,*) (b(ii,i),ii=1,13)
1538 write (iout,*) 'Type ',onelet(iloctyp(i))
1539 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1547 CCold(1,1,i)= b(7,i)
1548 CCold(2,2,i)=-b(7,i)
1549 CCold(2,1,i)= b(9,i)
1550 CCold(1,2,i)= b(9,i)
1551 CCold(1,1,-i)= b(7,i)
1552 CCold(2,2,-i)=-b(7,i)
1553 CCold(2,1,-i)=-b(9,i)
1554 CCold(1,2,-i)=-b(9,i)
1555 DDold(1,1,i)= b(6,i)
1556 DDold(2,2,i)=-b(6,i)
1557 DDold(2,1,i)= b(8,i)
1558 DDold(1,2,i)= b(8,i)
1559 DDold(1,1,-i)= b(6,i)
1560 DDold(2,2,-i)=-b(6,i)
1561 DDold(2,1,-i)=-b(8,i)
1562 DDold(1,2,-i)=-b(8,i)
1563 EEold(1,1,i)= b(10,i)+b(11,i)
1564 EEold(2,2,i)=-b(10,i)+b(11,i)
1565 EEold(2,1,i)= b(12,i)-b(13,i)
1566 EEold(1,2,i)= b(12,i)+b(13,i)
1567 EEold(1,1,-i)= b(10,i)+b(11,i)
1568 EEold(2,2,-i)=-b(10,i)+b(11,i)
1569 EEold(2,1,-i)=-b(12,i)+b(13,i)
1570 EEold(1,2,-i)=-b(12,i)-b(13,i)
1571 write(iout,*) "TU DOCHODZE"
1577 "Coefficients of the cumulants (independent of valence angles)"
1578 do i=-nloctyp+1,nloctyp-1
1579 write (iout,*) 'Type ',onelet(iloctyp(i))
1581 write(iout,'(2f10.5)') B(3,i),B(5,i)
1583 write(iout,'(2f10.5)') B(2,i),B(4,i)
1586 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1590 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1594 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1601 ! Read torsional parameters in old format
1603 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1605 read (itorp,*) ntortyp,nterm_old
1606 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1607 read (itorp,*) (itortyp(i),i=1,ntyp)
1609 !el from energy module--------
1610 allocate(v1(nterm_old,ntortyp,ntortyp))
1611 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1612 !el---------------------------
1618 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1624 write (iout,'(/a/)') 'Torsional constants:'
1627 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1628 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1636 ! Read torsional parameters
1638 IF (TOR_MODE.eq.0) THEN
1640 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1642 read (itorp,*) ntortyp
1643 read (itorp,*) (itortyp(i),i=1,ntyp)
1644 write (iout,*) 'ntortyp',ntortyp
1646 !el from energy module---------
1647 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1648 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1650 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1651 allocate(vlor2(maxlor,ntortyp,ntortyp))
1652 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1653 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1655 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1656 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1657 !el---------------------------
1659 do i=-ntortyp,ntortyp
1660 do j=-ntortyp,ntortyp
1666 !el---------------------------
1670 itortyp(i)=-itortyp(-i)
1672 ! write (iout,*) 'ntortyp',ntortyp
1674 do j=-ntortyp+1,ntortyp-1
1675 read (itorp,*) nterm(i,j,iblock),&
1677 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1678 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1681 do k=1,nterm(i,j,iblock)
1682 read (itorp,*) kk,v1(k,i,j,iblock),&
1684 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1685 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1686 v0ij=v0ij+si*v1(k,i,j,iblock)
1689 do k=1,nlor(i,j,iblock)
1690 read (itorp,*) kk,vlor1(k,i,j),&
1691 vlor2(k,i,j),vlor3(k,i,j)
1692 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1695 v0(-i,-j,iblock)=v0ij
1702 write (iout,'(/a/)') 'Torsional constants:'
1705 write (iout,*) 'ityp',i,' jtyp',j
1706 write (iout,*) 'Fourier constants'
1707 do k=1,nterm(i,j,iblock)
1708 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1711 write (iout,*) 'Lorenz constants'
1712 do k=1,nlor(i,j,iblock)
1713 write (iout,'(3(1pe15.5))') &
1714 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1721 ! 6/23/01 Read parameters for double torsionals
1723 !el from energy module------------
1724 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1725 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1726 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1727 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1728 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1729 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1730 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1731 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1732 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1733 !---------------------------------
1737 do j=-ntortyp+1,ntortyp-1
1738 do k=-ntortyp+1,ntortyp-1
1739 read (itordp,'(3a1)') t1,t2,t3
1740 ! write (iout,*) "OK onelett",
1743 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1744 .or. t3.ne.toronelet(k)) then
1745 write (iout,*) "Error in double torsional parameter file",&
1748 call MPI_Finalize(Ierror)
1750 stop "Error in double torsional parameter file"
1752 read (itordp,*) ntermd_1(i,j,k,iblock),&
1753 ntermd_2(i,j,k,iblock)
1754 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1755 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1756 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1757 ntermd_1(i,j,k,iblock))
1758 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1759 ntermd_1(i,j,k,iblock))
1760 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1761 ntermd_1(i,j,k,iblock))
1762 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1763 ntermd_1(i,j,k,iblock))
1764 ! Martix of D parameters for one dimesional foureir series
1765 do l=1,ntermd_1(i,j,k,iblock)
1766 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1767 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1768 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1769 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1770 ! write(iout,*) "whcodze" ,
1771 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1773 read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1774 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1775 v2s(m,l,i,j,k,iblock),&
1776 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1777 ! Martix of D parameters for two dimesional fourier series
1778 do l=1,ntermd_2(i,j,k,iblock)
1780 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1781 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1782 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1783 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1792 write (iout,*) 'Constants for double torsionals'
1795 do j=-ntortyp+1,ntortyp-1
1796 do k=-ntortyp+1,ntortyp-1
1797 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1798 ' nsingle',ntermd_1(i,j,k,iblock),&
1799 ' ndouble',ntermd_2(i,j,k,iblock)
1801 write (iout,*) 'Single angles:'
1802 do l=1,ntermd_1(i,j,k,iblock)
1803 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1804 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1805 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1806 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1809 write (iout,*) 'Pairs of angles:'
1810 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1811 do l=1,ntermd_2(i,j,k,iblock)
1812 write (iout,'(i5,20f10.5)') &
1813 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1816 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1817 do l=1,ntermd_2(i,j,k,iblock)
1818 write (iout,'(i5,20f10.5)') &
1819 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1820 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1830 itype2loc(i)=itortyp(i)
1833 ELSE IF (TOR_MODE.eq.1) THEN
1835 !C read valence-torsional parameters
1836 read (itorp,*) ntortyp
1838 write (iout,*) "Valence-torsional parameters read in ntortyp",&
1840 read (itorp,*) (itortyp(i),i=1,ntyp)
1841 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1844 itype2loc(i)=itortyp(i)
1848 itortyp(i)=-itortyp(-i)
1850 do i=-ntortyp+1,ntortyp-1
1851 do j=-ntortyp+1,ntortyp-1
1852 !C first we read the cos and sin gamma parameters
1853 read (itorp,'(13x,a)') string
1854 write (iout,*) i,j,string
1856 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1857 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1858 do k=1,nterm_kcc(j,i)
1859 do l=1,nterm_kcc_Tb(j,i)
1860 do ll=1,nterm_kcc_Tb(j,i)
1861 read (itorp,*) ii,jj,kk, &
1862 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1870 !c AL 4/8/16: Calculate coefficients from one-body parameters
1872 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1873 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
1874 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
1875 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1876 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1879 print *,i,itortyp(i)
1880 itortyp(i)=itype2loc(i)
1883 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
1885 do i=-ntortyp+1,ntortyp-1
1886 do j=-ntortyp+1,ntortyp-1
1889 do k=1,nterm_kcc_Tb(j,i)
1890 do l=1,nterm_kcc_Tb(j,i)
1891 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
1892 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1893 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
1894 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1897 do k=1,nterm_kcc_Tb(j,i)
1898 do l=1,nterm_kcc_Tb(j,i)
1900 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1901 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1902 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1903 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1905 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1906 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1907 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1908 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1912 !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)
1916 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1920 if (tor_mode.gt.0 .and. lprint) then
1921 !c Print valence-torsional parameters
1922 write (iout,'(a)') &
1923 "Parameters of the valence-torsional potentials"
1924 do i=-ntortyp+1,ntortyp-1
1925 do j=-ntortyp+1,ntortyp-1
1926 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1927 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1928 do k=1,nterm_kcc(j,i)
1929 do l=1,nterm_kcc_Tb(j,i)
1930 do ll=1,nterm_kcc_Tb(j,i)
1931 write (iout,'(3i5,2f15.4)')&
1932 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1941 !elwrite(iout,*) "parmread kontrol sc-bb"
1942 ! Read of Side-chain backbone correlation parameters
1943 ! Modified 11 May 2012 by Adasko
1946 read (isccor,*) nsccortyp
1949 !c maxinter is maximum interaction sites
1950 !write(iout,*)"maxterm_sccor",maxterm_sccor
1951 !el from module energy-------------
1952 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1953 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1954 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1955 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1956 !-----------------------------------
1957 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1958 !-----------------------------------
1959 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1960 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1961 -nsccortyp:nsccortyp))
1962 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1963 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1964 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1965 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1966 !-----------------------------------
1967 do i=-nsccortyp,nsccortyp
1968 do j=-nsccortyp,nsccortyp
1972 !-----------------------------------
1974 read (isccor,*) (isccortyp(i),i=1,ntyp)
1976 isccortyp(i)=-isccortyp(-i)
1978 iscprol=isccortyp(20)
1979 ! write (iout,*) 'ntortyp',ntortyp
1981 !c maxinter is maximum interaction sites
1986 nterm_sccor(i,j),nlor_sccor(i,j)
1992 nterm_sccor(-i,j)=nterm_sccor(i,j)
1993 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1994 nterm_sccor(i,-j)=nterm_sccor(i,j)
1995 do k=1,nterm_sccor(i,j)
1996 read (isccor,*) kk,v1sccor(k,l,i,j),&
1998 if (j.eq.iscprol) then
1999 if (i.eq.isccortyp(10)) then
2000 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2001 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2003 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2004 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2005 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2006 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2007 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2008 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2009 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2010 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2013 if (i.eq.isccortyp(10)) then
2014 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2015 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2017 if (j.eq.isccortyp(10)) then
2018 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2019 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2021 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2022 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2023 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2024 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2025 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2026 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2030 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2031 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2032 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2033 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2036 do k=1,nlor_sccor(i,j)
2037 read (isccor,*) kk,vlor1sccor(k,i,j),&
2038 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2039 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2040 (1+vlor3sccor(k,i,j)**2)
2042 v0sccor(l,i,j)=v0ijsccor
2043 v0sccor(l,-i,j)=v0ijsccor1
2044 v0sccor(l,i,-j)=v0ijsccor2
2045 v0sccor(l,-i,-j)=v0ijsccor3
2051 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
2054 write (iout,*) 'ityp',i,' jtyp',j
2055 write (iout,*) 'Fourier constants'
2056 do k=1,nterm_sccor(i,j)
2057 write (iout,'(2(1pe15.5))') &
2058 (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
2060 write (iout,*) 'Lorenz constants'
2061 do k=1,nlor_sccor(i,j)
2062 write (iout,'(3(1pe15.5))') &
2063 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2070 ! Read electrostatic-interaction parameters
2073 write (iout,'(/a)') 'Electrostatic interaction constants:'
2074 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2075 'IT','JT','APP','BPP','AEL6','AEL3'
2077 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
2078 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
2079 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
2080 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
2085 app (i,j)=epp(i,j)*rri*rri
2086 bpp (i,j)=-2.0D0*epp(i,j)*rri
2087 ael6(i,j)=elpp6(i,j)*4.2D0**6
2088 ael3(i,j)=elpp3(i,j)*4.2D0**3
2089 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2094 ! Read side-chain interaction parameters.
2096 !el from module energy - COMMON.INTERACT-------
2097 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2098 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2099 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2100 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2101 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2102 allocate(epslip(ntyp,ntyp))
2113 !--------------------------------
2115 read (isidep,*) ipot,expon
2116 !el if (ipot.lt.1 .or. ipot.gt.5) then
2117 ! write (iout,'(2a)') 'Error while reading SC interaction',&
2118 ! 'potential file - unknown potential type.'
2122 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2123 ', exponents are ',expon,2*expon
2124 ! goto (10,20,30,30,40) ipot
2126 !----------------------- LJ potential ---------------------------------
2128 ! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2129 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2131 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2132 write (iout,'(a/)') 'The epsilon array:'
2133 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2134 write (iout,'(/a)') 'One-body parameters:'
2135 write (iout,'(a,4x,a)') 'residue','sigma'
2136 write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
2139 !----------------------- LJK potential --------------------------------
2141 ! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2142 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2143 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2145 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2146 write (iout,'(a/)') 'The epsilon array:'
2147 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2148 write (iout,'(/a)') 'One-body parameters:'
2149 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2150 write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2154 !---------------------- GB or BP potential -----------------------------
2157 if (scelemode.eq.0) then
2159 read (isidep,*)(eps(i,j),j=i,ntyp)
2161 read (isidep,*)(sigma0(i),i=1,ntyp)
2162 read (isidep,*)(sigii(i),i=1,ntyp)
2163 read (isidep,*)(chip(i),i=1,ntyp)
2164 read (isidep,*)(alp(i),i=1,ntyp)
2166 read (isidep,*)(epslip(i,j),j=i,ntyp)
2168 ! For the GB potential convert sigma'**2 into chi'
2171 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2175 write (iout,'(/a/)') 'Parameters of the BP potential:'
2176 write (iout,'(a/)') 'The epsilon array:'
2177 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2178 write (iout,'(/a)') 'One-body parameters:'
2179 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2181 write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
2182 chip(i),alp(i),i=1,ntyp)
2185 allocate(icharge(ntyp1))
2186 ! print *,ntyp,icharge(i)
2188 read (isidep,*) (icharge(i),i=1,ntyp)
2189 print *,ntyp,icharge(i)
2190 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2191 write (2,*) "icharge",(icharge(i),i=1,ntyp)
2192 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2193 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2194 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2195 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2196 allocate(epsintab(ntyp,ntyp))
2197 allocate(dtail(2,ntyp,ntyp))
2198 print *,"control line 1"
2199 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2200 allocate(wqdip(2,ntyp,ntyp))
2201 allocate(wstate(4,ntyp,ntyp))
2202 allocate(dhead(2,2,ntyp,ntyp))
2203 allocate(nstate(ntyp,ntyp))
2204 allocate(debaykap(ntyp,ntyp))
2205 print *,"control line 2"
2206 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2207 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2211 ! write (*,*) "Im in ALAB", i, " ", j
2213 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2214 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2215 chis(i,j),chis(j,i), &
2216 nstate(i,j),(wstate(k,i,j),k=1,4), &
2217 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2218 dtail(1,i,j),dtail(2,i,j), &
2219 epshead(i,j),sig0head(i,j), &
2220 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2221 alphapol(i,j),alphapol(j,i), &
2222 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
2223 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2229 sigma(i,j) = sigma(j,i)
2230 sigmap1(i,j)=sigmap1(j,i)
2231 sigmap2(i,j)=sigmap2(j,i)
2232 sigiso1(i,j)=sigiso1(j,i)
2233 sigiso2(i,j)=sigiso2(j,i)
2234 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2235 nstate(i,j) = nstate(j,i)
2236 dtail(1,i,j) = dtail(1,j,i)
2237 dtail(2,i,j) = dtail(2,j,i)
2239 alphasur(k,i,j) = alphasur(k,j,i)
2240 wstate(k,i,j) = wstate(k,j,i)
2241 alphiso(k,i,j) = alphiso(k,j,i)
2244 dhead(2,1,i,j) = dhead(1,1,j,i)
2245 dhead(2,2,i,j) = dhead(1,2,j,i)
2246 dhead(1,1,i,j) = dhead(2,1,j,i)
2247 dhead(1,2,i,j) = dhead(2,2,j,i)
2249 epshead(i,j) = epshead(j,i)
2250 sig0head(i,j) = sig0head(j,i)
2253 wqdip(k,i,j) = wqdip(k,j,i)
2256 wquad(i,j) = wquad(j,i)
2257 epsintab(i,j) = epsintab(j,i)
2258 debaykap(i,j)=debaykap(j,i)
2259 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2265 !--------------------- GBV potential -----------------------------------
2267 ! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2268 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2269 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2270 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2272 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2273 write (iout,'(a/)') 'The epsilon array:'
2274 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2275 write (iout,'(/a)') 'One-body parameters:'
2276 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2277 's||/s_|_^2',' chip ',' alph '
2278 write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2279 sigii(i),chip(i),alp(i),i=1,ntyp)
2282 write (iout,'(2a)') 'Error while reading SC interaction',&
2283 'potential file - unknown potential type.'
2289 !-----------------------------------------------------------------------
2290 ! Calculate the "working" parameters of SC interactions.
2292 !el from module energy - COMMON.INTERACT-------
2293 ! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2294 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1))
2295 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) !(ntyp,ntyp)
2296 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2297 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2298 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2299 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2301 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2309 if (scelemode.eq.0) then
2316 !--------------------------------
2321 epslip(i,j)=epslip(j,i)
2324 if (scelemode.eq.0) then
2327 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2328 sigma(j,i)=sigma(i,j)
2329 rs0(i,j)=dwa16*sigma(i,j)
2334 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2335 'Working parameters of the SC interactions:',&
2336 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2341 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2351 sigeps=dsign(1.0D0,epsij)
2353 aa_aq(i,j)=epsij*rrij*rrij
2354 bb_aq(i,j)=-sigeps*epsij*rrij
2355 aa_aq(j,i)=aa_aq(i,j)
2356 bb_aq(j,i)=bb_aq(i,j)
2357 epsijlip=epslip(i,j)
2358 sigeps=dsign(1.0D0,epsijlip)
2359 epsijlip=dabs(epsijlip)
2360 aa_lip(i,j)=epsijlip*rrij*rrij
2361 bb_lip(i,j)=-sigeps*epsijlip*rrij
2362 aa_lip(j,i)=aa_lip(i,j)
2363 bb_lip(j,i)=bb_lip(i,j)
2364 if ((ipot.gt.2).and. (scelemode.eq.0))then
2365 sigt1sq=sigma0(i)**2
2366 sigt2sq=sigma0(j)**2
2369 ratsig1=sigt2sq/sigt1sq
2370 ratsig2=1.0D0/ratsig1
2371 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2372 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2373 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2377 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2378 sigmaii(i,j)=rsum_max
2379 sigmaii(j,i)=rsum_max
2381 ! sigmaii(i,j)=r0(i,j)
2382 ! sigmaii(j,i)=r0(i,j)
2384 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2385 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2386 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2387 augm(i,j)=epsij*r_augm**(2*expon)
2388 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2395 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2396 restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2397 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2401 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2402 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2403 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2404 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2405 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2406 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2407 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2408 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2409 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2410 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2411 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2412 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2413 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2414 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2415 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2416 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2425 read (isidep_nucl,*) ipot_nucl
2426 ! print *,"TU?!",ipot_nucl
2427 if (ipot_nucl.eq.1) then
2428 do i=1,ntyp_molec(2)
2429 do j=i,ntyp_molec(2)
2430 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2431 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2435 do i=1,ntyp_molec(2)
2436 do j=i,ntyp_molec(2)
2437 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2438 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2439 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2443 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2444 do i=1,ntyp_molec(2)
2445 do j=i,ntyp_molec(2)
2446 rrij=sigma_nucl(i,j)
2450 epsij=4*eps_nucl(i,j)
2451 sigeps=dsign(1.0D0,epsij)
2453 aa_nucl(i,j)=epsij*rrij*rrij
2454 bb_nucl(i,j)=-sigeps*epsij*rrij
2455 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2456 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2457 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2458 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2459 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2460 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2463 aa_nucl(i,j)=aa_nucl(j,i)
2464 bb_nucl(i,j)=bb_nucl(j,i)
2465 ael3_nucl(i,j)=ael3_nucl(j,i)
2466 ael6_nucl(i,j)=ael6_nucl(j,i)
2467 ael63_nucl(i,j)=ael63_nucl(j,i)
2468 ael32_nucl(i,j)=ael32_nucl(j,i)
2469 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2470 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2471 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2472 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2473 eps_nucl(i,j)=eps_nucl(j,i)
2474 sigma_nucl(i,j)=sigma_nucl(j,i)
2475 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2479 write(iout,*) "tube param"
2480 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2481 ccavtubpep,dcavtubpep,tubetranenepep
2482 sigmapeptube=sigmapeptube**6
2483 sigeps=dsign(1.0D0,epspeptube)
2484 epspeptube=dabs(epspeptube)
2485 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2486 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2487 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2489 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2490 ccavtub(i),dcavtub(i),tubetranene(i)
2491 sigmasctube=sigmasctube**6
2492 sigeps=dsign(1.0D0,epssctube)
2493 epssctube=dabs(epssctube)
2494 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2495 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2496 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2498 !-----------------READING SC BASE POTENTIALS-----------------------------
2499 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2500 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2501 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2502 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2503 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2504 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2505 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2506 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2507 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2508 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2509 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2510 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2511 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2512 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2513 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2514 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2517 do i=1,ntyp_molec(1)
2518 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2519 write (*,*) "Im in ", i, " ", j
2520 read(isidep_scbase,*) &
2521 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2522 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2523 write(*,*) "eps",eps_scbase(i,j)
2524 read(isidep_scbase,*) &
2525 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2526 chis_scbase(i,j,1),chis_scbase(i,j,2)
2527 read(isidep_scbase,*) &
2528 dhead_scbasei(i,j), &
2529 dhead_scbasej(i,j), &
2530 rborn_scbasei(i,j),rborn_scbasej(i,j)
2531 read(isidep_scbase,*) &
2532 (wdipdip_scbase(k,i,j),k=1,3), &
2533 (wqdip_scbase(k,i,j),k=1,2)
2534 read(isidep_scbase,*) &
2535 alphapol_scbase(i,j), &
2536 epsintab_scbase(i,j)
2539 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2540 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2542 do i=1,ntyp_molec(1)
2543 do j=1,ntyp_molec(2)-1
2544 epsij=eps_scbase(i,j)
2545 rrij=sigma_scbase(i,j)
2550 sigeps=dsign(1.0D0,epsij)
2552 aa_scbase(i,j)=epsij*rrij*rrij
2553 bb_scbase(i,j)=-sigeps*epsij*rrij
2556 !-----------------READING PEP BASE POTENTIALS-------------------
2557 allocate(eps_pepbase(ntyp_molec(2)))
2558 allocate(sigma_pepbase(ntyp_molec(2)))
2559 allocate(chi_pepbase(ntyp_molec(2),2))
2560 allocate(chipp_pepbase(ntyp_molec(2),2))
2561 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2562 allocate(sigmap1_pepbase(ntyp_molec(2)))
2563 allocate(sigmap2_pepbase(ntyp_molec(2)))
2564 allocate(chis_pepbase(ntyp_molec(2),2))
2565 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2568 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2569 write (*,*) "Im in ", i, " ", j
2570 read(isidep_pepbase,*) &
2571 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2572 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2573 write(*,*) "eps",eps_pepbase(j)
2574 read(isidep_pepbase,*) &
2575 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2576 chis_pepbase(j,1),chis_pepbase(j,2)
2577 read(isidep_pepbase,*) &
2578 (wdipdip_pepbase(k,j),k=1,3)
2580 allocate(aa_pepbase(ntyp_molec(2)))
2581 allocate(bb_pepbase(ntyp_molec(2)))
2583 do j=1,ntyp_molec(2)-1
2584 epsij=eps_pepbase(j)
2585 rrij=sigma_pepbase(j)
2590 sigeps=dsign(1.0D0,epsij)
2592 aa_pepbase(j)=epsij*rrij*rrij
2593 bb_pepbase(j)=-sigeps*epsij*rrij
2595 !--------------READING SC PHOSPHATE-------------------------------------
2596 !--------------READING SC PHOSPHATE-------------------------------------
2597 allocate(eps_scpho(ntyp_molec(1)))
2598 allocate(sigma_scpho(ntyp_molec(1)))
2599 allocate(chi_scpho(ntyp_molec(1),2))
2600 allocate(chipp_scpho(ntyp_molec(1),2))
2601 allocate(alphasur_scpho(4,ntyp_molec(1)))
2602 allocate(sigmap1_scpho(ntyp_molec(1)))
2603 allocate(sigmap2_scpho(ntyp_molec(1)))
2604 allocate(chis_scpho(ntyp_molec(1),2))
2605 allocate(wqq_scpho(ntyp_molec(1)))
2606 allocate(wqdip_scpho(2,ntyp_molec(1)))
2607 allocate(alphapol_scpho(ntyp_molec(1)))
2608 allocate(epsintab_scpho(ntyp_molec(1)))
2609 allocate(dhead_scphoi(ntyp_molec(1)))
2610 allocate(rborn_scphoi(ntyp_molec(1)))
2611 allocate(rborn_scphoj(ntyp_molec(1)))
2612 allocate(alphi_scpho(ntyp_molec(1)))
2616 do j=1,ntyp_molec(1) ! without U then we will take T for U
2617 write (*,*) "Im in scpho ", i, " ", j
2618 read(isidep_scpho,*) &
2619 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2620 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2621 write(*,*) "eps",eps_scpho(j)
2622 read(isidep_scpho,*) &
2623 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2624 chis_scpho(j,1),chis_scpho(j,2)
2625 read(isidep_scpho,*) &
2626 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2627 read(isidep_scpho,*) &
2628 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2632 allocate(aa_scpho(ntyp_molec(1)))
2633 allocate(bb_scpho(ntyp_molec(1)))
2634 do j=1,ntyp_molec(1)
2641 sigeps=dsign(1.0D0,epsij)
2643 aa_scpho(j)=epsij*rrij*rrij
2644 bb_scpho(j)=-sigeps*epsij*rrij
2648 read(isidep_peppho,*) &
2649 eps_peppho,sigma_peppho
2650 read(isidep_peppho,*) &
2651 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2652 read(isidep_peppho,*) &
2653 (wqdip_peppho(k),k=1,2)
2661 sigeps=dsign(1.0D0,epsij)
2663 aa_peppho=epsij*rrij*rrij
2664 bb_peppho=-sigeps*epsij*rrij
2667 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2675 ! Define the SC-p interaction constants
2684 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2685 read (itorp_nucl,*) ntortyp_nucl
2686 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2687 !el from energy module---------
2688 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2689 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2691 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2692 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2693 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2694 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2696 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2697 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2698 !el---------------------------
2701 !el--------------------
2702 read (itorp_nucl,*) &
2703 (itortyp_nucl(i),i=1,ntyp_molec(2))
2704 ! print *,itortyp_nucl(:)
2705 !c write (iout,*) 'ntortyp',ntortyp
2708 read (itorp_nucl,*) nterm_nucl(i,j),nlor_nucl(i,j)
2709 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2712 do k=1,nterm_nucl(i,j)
2713 read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2714 v0ij=v0ij+si*v1_nucl(k,i,j)
2717 do k=1,nlor_nucl(i,j)
2718 read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),&
2719 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2720 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2727 !elwrite(iout,*) "parmread kontrol before oldscp"
2729 ! Define the SC-p interaction constants
2733 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2735 ! aad(i,1)=0.3D0*4.0D0**12
2736 ! Following line for constants currently implemented
2737 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2738 aad(i,1)=1.5D0*4.0D0**12
2739 ! aad(i,1)=0.17D0*5.6D0**12
2741 ! "Soft" SC-p repulsion
2743 ! Following line for constants currently implemented
2744 ! aad(i,1)=0.3D0*4.0D0**6
2745 ! "Hard" SC-p repulsion
2746 bad(i,1)=3.0D0*4.0D0**6
2747 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2756 ! 8/9/01 Read the SC-p interaction constants from file
2759 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
2762 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2763 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2764 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2765 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2769 write (iout,*) "Parameters of SC-p interactions:"
2771 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2772 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2776 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2778 do i=1,ntyp_molec(2)
2779 read (iscpp_nucl,*) eps_scp_nucl(i),rscp_nucl(i)
2781 do i=1,ntyp_molec(2)
2782 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2783 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2785 r0pp=1.12246204830937298142*5.16158
2791 ! Define the constants of the disulfide bridge
2795 ! Old arbitrary potential - commented out.
2800 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2801 ! energy surface of diethyl disulfide.
2802 ! A. Liwo and U. Kozlowska, 11/24/03
2813 write (iout,'(/a)') "Disulfide bridge parameters:"
2814 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2815 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2816 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2817 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2821 end subroutine parmread
2823 !-----------------------------------------------------------------------------
2825 !-----------------------------------------------------------------------------
2826 subroutine mygetenv(string,var)
2830 ! This subroutine passes the environmental variables to FORTRAN program.
2831 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
2832 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
2833 ! reads the environmental variables from $HOME/.env
2835 ! Usage: As for the standard FORTRAN GETENV subroutine.
2837 ! Purpose: some versions/installations of MPI do not transfer the environmental
2838 ! variables to slave processors, if these variables are set in the shell script
2839 ! from which mpirun is called.
2848 character*(*) :: string,var
2849 #if defined(MYGETENV) && defined(MPI)
2850 ! include "DIMENSIONS.ZSCOPT"
2852 ! include "COMMON.MPI"
2853 !el character*360 ucase
2855 character(len=360) :: string1(360),karta
2856 character(len=240) :: home
2859 call getenv("HOME",home)
2860 open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
2862 read (99,end=111,err=111,'(a)') karta
2866 call split_string(karta,string1,80,n)
2867 if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
2868 string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
2870 print *,"Processor",me,": ",var(:ilen(var)),&
2871 " assigned to ",string(:ilen(string))
2876 111 print *,"Environment variable ",string(:ilen(string))," not set."
2879 112 print *,"Error opening environment file!"
2881 call getenv(string,var)
2884 end subroutine mygetenv
2885 !-----------------------------------------------------------------------------
2887 !-----------------------------------------------------------------------------
2888 subroutine read_general_data(*)
2890 use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
2891 scelemode,TUBEmode,tor_mode
2893 use energy_data, only:distchainmax,tubeR0,tubecenter
2894 use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
2895 bordtubebot,tubebufthick,buftubebot,buftubetop
2897 ! include "DIMENSIONS"
2898 ! include "DIMENSIONS.ZSCOPT"
2899 ! include "DIMENSIONS.FREE"
2900 ! include "COMMON.TORSION"
2901 ! include "COMMON.INTERACT"
2902 ! include "COMMON.IOUNITS"
2903 ! include "COMMON.TIME1"
2904 ! include "COMMON.PROT"
2905 ! include "COMMON.PROTFILES"
2906 ! include "COMMON.CHAIN"
2907 ! include "COMMON.NAMES"
2908 ! include "COMMON.FFIELD"
2909 ! include "COMMON.ENEPS"
2910 ! include "COMMON.WEIGHTS"
2911 ! include "COMMON.FREE"
2912 ! include "COMMON.CONTROL"
2913 ! include "COMMON.ENERGIES"
2914 character(len=800) :: controlcard
2915 integer :: i,j,k,ii,n_ene_found
2916 integer :: ind,itype1,itype2,itypf,itypsc,itypp
2919 !el character*16 ucase
2920 character(len=16) :: key
2922 call card_concat(controlcard,.true.)
2923 call readi(controlcard,"N_ENE",n_eneW,max_eneW)
2924 if (n_eneW.gt.max_eneW) then
2925 write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
2929 call readi(controlcard,"NPARMSET",nparmset,1)
2930 !elwrite(iout,*)"in read_gen data"
2931 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
2932 call readi(controlcard,"IPARMPRINT",iparmprint,1)
2933 write (iout,*) "PARMPRINT",iparmprint
2934 if (nparmset.gt.max_parm) then
2935 write (iout,*) "Error: parameter out of range: NPARMSET",&
2939 !elwrite(iout,*)"in read_gen data"
2940 call readi(controlcard,"MAXIT",maxit,5000)
2941 call reada(controlcard,"FIMIN",fimin,1.0d-3)
2942 call readi(controlcard,"ENSEMBLES",ensembles,0)
2943 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
2944 write (iout,*) "Number of energy parameter sets",nparmset
2945 allocate(isampl(nparmset))
2946 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
2947 write (iout,*) "MaxSlice",MaxSlice
2948 call readi(controlcard,"NSLICE",nslice,1)
2949 !elwrite(iout,*)"in read_gen data"
2951 if (nslice.gt.MaxSlice) then
2952 write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
2956 write (iout,*) "Frequency of storing conformations",&
2957 (isampl(i),i=1,nparmset)
2958 write (iout,*) "Maxit",maxit," Fimin",fimin
2959 call readi(controlcard,"NQ",nQ,1)
2960 if (nQ.gt.MaxQ) then
2961 write (iout,*) "Error: parameter out of range: NQ",nq,&
2966 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
2967 call reada(controlcard,"DELTA",delta,1.0d-2)
2968 call readi(controlcard,"EINICHECK",einicheck,2)
2969 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
2970 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
2971 call readi(controlcard,"RESCALE",rescale_modeW,1)
2972 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2973 call reada(controlcard,'BOXY',boxysize,100.0d0)
2974 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2975 call readi(controlcard,"SCELEMODE",scelemode,0)
2976 print *,"SCELE",scelemode
2977 call readi(controlcard,'TORMODE',tor_mode,0)
2978 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2979 write(iout,*) "torsional and valence angle mode",tor_mode
2981 call readi(controlcard,'TUBEMOD',tubemode,0)
2983 if (TUBEmode.gt.0) then
2984 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
2985 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
2986 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
2987 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
2988 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
2989 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
2990 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
2991 buftubebot=bordtubebot+tubebufthick
2992 buftubetop=bordtubetop-tubebufthick
2994 ions=index(controlcard,"IONS").gt.0
2995 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2996 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2997 write(iout,*) "R_CUT_ELE=",r_cut_ele
2998 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
2999 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
3000 call readi(controlcard,'SYM',symetr,1)
3001 write (iout,*) "DISTCHAINMAX",distchainmax
3002 write (iout,*) "delta",delta
3003 write (iout,*) "einicheck",einicheck
3004 write (iout,*) "rescale_mode",rescale_modeW
3006 bxfile=index(controlcard,"BXFILE").gt.0
3007 cxfile=index(controlcard,"CXFILE").gt.0
3008 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
3010 histfile=index(controlcard,"HISTFILE").gt.0
3011 histout=index(controlcard,"HISTOUT").gt.0
3012 entfile=index(controlcard,"ENTFILE").gt.0
3013 zscfile=index(controlcard,"ZSCFILE").gt.0
3014 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
3015 write (iout,*) "with_dihed_constr ",with_dihed_constr
3016 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3018 end subroutine read_general_data
3019 !------------------------------------------------------------------------------
3020 subroutine read_efree(*)
3022 ! Read molecular data
3025 ! include 'DIMENSIONS'
3026 ! include 'DIMENSIONS.ZSCOPT'
3027 ! include 'DIMENSIONS.COMPAR'
3028 ! include 'DIMENSIONS.FREE'
3029 ! include 'COMMON.IOUNITS'
3030 ! include 'COMMON.TIME1'
3031 ! include 'COMMON.SBRIDGE'
3032 ! include 'COMMON.CONTROL'
3033 ! include 'COMMON.CHAIN'
3034 ! include 'COMMON.HEADER'
3035 ! include 'COMMON.GEO'
3036 ! include 'COMMON.FREE'
3037 character(len=320) :: controlcard !,ucase
3038 integer :: iparm,ib,i,j,npars
3048 ! call alloc_wham_arrays
3049 ! allocate(nT_h(nParmSet))
3050 ! allocate(replica(nParmSet))
3051 ! allocate(umbrella(nParmSet))
3052 ! allocate(read_iset(nParmSet))
3053 ! allocate(nT_h(nParmSet))
3057 call card_concat(controlcard,.true.)
3058 call readi(controlcard,'NT',nT_h(iparm),1)
3059 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
3061 if (nT_h(iparm).gt.MaxT_h) then
3062 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
3066 replica(iparm)=index(controlcard,"REPLICA").gt.0
3067 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
3068 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
3069 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
3070 replica(iparm)," umbrella ",umbrella(iparm),&
3071 " read_iset",read_iset(iparm)
3074 call card_concat(controlcard,.true.)
3075 call readi(controlcard,'NR',nR(ib,iparm),1)
3076 if (umbrella(iparm)) then
3079 nRR(ib,iparm)=nR(ib,iparm)
3081 if (nR(ib,iparm).gt.MaxR) then
3082 write (iout,*) "Error: parameter out of range: NR",&
3086 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
3087 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
3088 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
3091 call card_concat(controlcard,.true.)
3092 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
3094 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
3099 write (iout,*) "ib",ib," beta_h",&
3100 1.0d0/(0.001987*beta_h(ib,iparm))
3101 write (iout,*) "nR",nR(ib,iparm)
3102 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
3104 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
3105 "q0",(q0(j,i,ib,iparm),j=1,nQ)
3117 nR(ib,iparm)=nR(ib,1)
3118 if (umbrella(iparm)) then
3121 nRR(ib,iparm)=nR(ib,1)
3123 beta_h(ib,iparm)=beta_h(ib,1)
3125 f(i,ib,iparm)=f(i,ib,1)
3127 KH(j,i,ib,iparm)=KH(j,i,ib,1)
3128 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
3131 replica(iparm)=replica(1)
3132 umbrella(iparm)=umbrella(1)
3133 read_iset(iparm)=read_iset(1)
3140 end subroutine read_efree
3141 !-----------------------------------------------------------------------------
3142 subroutine read_protein_data(*)
3144 ! include "DIMENSIONS"
3145 ! include "DIMENSIONS.ZSCOPT"
3146 ! include "DIMENSIONS.FREE"
3150 integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
3151 ! include "COMMON.MPI"
3153 ! include "COMMON.CHAIN"
3154 ! include "COMMON.IOUNITS"
3155 ! include "COMMON.PROT"
3156 ! include "COMMON.PROTFILES"
3157 ! include "COMMON.NAMES"
3158 ! include "COMMON.FREE"
3159 ! include "COMMON.OBCINKA"
3160 character(len=64) :: nazwa
3161 character(len=16000) :: controlcard
3162 integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
3163 !el external ilen,iroof
3172 ! Read names of files with conformation data.
3173 if (replica(iparm)) then
3179 do ii=1,nRR(ib,iparm)
3180 write (iout,*) "Parameter set",iparm," temperature",ib,&
3183 call card_concat(controlcard,.true.)
3184 write (iout,*) controlcard(:ilen(controlcard))
3185 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
3186 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
3187 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
3188 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
3189 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
3190 maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
3191 call reada(controlcard,"TIME_START",&
3192 time_start_collect(ii,ib,iparm),0.0d0)
3193 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
3195 write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
3196 " rec_end",rec_end(ii,ib,iparm)
3197 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
3198 " time_end",time_end_collect(ii,ib,iparm)
3200 if (replica(iparm)) then
3201 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
3202 write (iout,*) "Number of trajectories",totraj(ii,iparm)
3205 if (nfile_bin(ii,ib,iparm).lt.2 &
3206 .and. nfile_asc(ii,ib,iparm).eq.0 &
3207 .and. nfile_cx(ii,ib,iparm).eq.0) then
3208 write (iout,*) "Error - no action specified!"
3211 if (nfile_bin(ii,ib,iparm).gt.0) then
3212 call card_concat(controlcard,.false.)
3213 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
3214 maxfile_prot,nfile_bin(ii,ib,iparm))
3216 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
3217 write(iout,*) (protfiles(i,1,ii,ib,iparm),&
3218 i=1,nfile_bin(ii,ib,iparm))
3221 if (nfile_asc(ii,ib,iparm).gt.0) then
3222 call card_concat(controlcard,.false.)
3223 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3224 maxfile_prot,nfile_asc(ii,ib,iparm))
3226 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
3227 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3228 i=1,nfile_asc(ii,ib,iparm))
3230 else if (nfile_cx(ii,ib,iparm).gt.0) then
3231 call card_concat(controlcard,.false.)
3232 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3233 maxfile_prot,nfile_cx(ii,ib,iparm))
3235 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
3236 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3237 i=1,nfile_cx(ii,ib,iparm))
3246 end subroutine read_protein_data
3247 !-------------------------------------------------------------------------------
3248 subroutine readsss(rekord,lancuch,wartosc,default)
3250 character*(*) :: rekord,lancuch,wartosc,default
3251 character(len=80) :: aux
3252 integer :: lenlan,lenrec,iread,ireade
3256 lenlan=ilen(lancuch)
3258 iread=index(rekord,lancuch(:lenlan)//"=")
3259 ! print *,"rekord",rekord," lancuch",lancuch
3260 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3261 if (iread.eq.0) then
3265 iread=iread+lenlan+1
3266 ! print *,"iread",iread
3267 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3268 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3270 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3272 ! print *,"iread",iread
3273 if (iread.gt.lenrec) then
3278 ! print *,"ireade",ireade
3279 do while (ireade.lt.lenrec .and. &
3280 .not.iblnk(rekord(ireade:ireade)))
3283 wartosc=rekord(iread:ireade)
3285 end subroutine readsss
3286 !----------------------------------------------------------------------------
3287 subroutine multreads(rekord,lancuch,tablica,dim,default)
3290 character*(*) rekord,lancuch,tablica(dim),default
3291 character(len=80) :: aux
3292 integer :: lenlan,lenrec,iread,ireade
3299 lenlan=ilen(lancuch)
3301 iread=index(rekord,lancuch(:lenlan)//"=")
3302 ! print *,"rekord",rekord," lancuch",lancuch
3303 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3304 if (iread.eq.0) return
3305 iread=iread+lenlan+1
3307 ! print *,"iread",iread
3308 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3309 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3311 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3313 ! print *,"iread",iread
3314 if (iread.gt.lenrec) return
3316 ! print *,"ireade",ireade
3317 do while (ireade.lt.lenrec .and. &
3318 .not.iblnk(rekord(ireade:ireade)))
3321 tablica(i)=rekord(iread:ireade)
3324 end subroutine multreads
3325 !----------------------------------------------------------------------------
3326 subroutine split_string(rekord,tablica,dim,nsub)
3328 integer :: dim,nsub,i,ii,ll,kk
3329 character*(*) tablica(dim)
3330 character*(*) rekord
3340 ! Find the start of term name
3342 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
3345 ! Parse the name into TABLICA(i) until blank found
3346 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
3348 tablica(i)(kk:kk)=rekord(ii:ii)
3351 if (kk.gt.0) nsub=nsub+1
3352 if (ii.gt.ll) return
3355 end subroutine split_string
3356 !--------------------------------------------------------------------------------
3358 !--------------------------------------------------------------------------------
3359 subroutine read_compar
3361 ! Read molecular data
3363 use conform_compar, only:alloc_compar_arrays
3364 use control_data, only:pdbref
3365 use geometry_data, only:deg2rad,rad2deg
3367 ! include 'DIMENSIONS'
3368 ! include 'DIMENSIONS.ZSCOPT'
3369 ! include 'DIMENSIONS.COMPAR'
3370 ! include 'DIMENSIONS.FREE'
3371 ! include 'COMMON.IOUNITS'
3372 ! include 'COMMON.TIME1'
3373 ! include 'COMMON.SBRIDGE'
3374 ! include 'COMMON.CONTROL'
3375 ! include 'COMMON.COMPAR'
3376 ! include 'COMMON.CHAIN'
3377 ! include 'COMMON.HEADER'
3378 ! include 'COMMON.GEO'
3379 ! include 'COMMON.FREE'
3380 character(len=320) :: controlcard !,ucase
3381 character(len=64) :: wfile
3385 !elwrite(iout,*)"jestesmy w read_compar"
3386 call card_concat(controlcard,.true.)
3387 pdbref=(index(controlcard,'PDBREF').gt.0)
3388 call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
3389 call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
3390 call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
3391 call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
3392 verbose = index(controlcard,"VERBOSE").gt.0
3393 lgrp=index(controlcard,"STATIN").gt.0
3394 lgrp_out=index(controlcard,"STATOUT").gt.0
3395 merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
3396 binary = index(controlcard,"BINARY").gt.0
3397 rmscut_base_up=rmscut_base_up/50
3398 rmscut_base_low=rmscut_base_low/50
3399 call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
3400 call readi(controlcard,'NLEVEL',nlevel,1)
3401 if (nlevel.lt.0) then
3403 call alloc_compar_arrays(maxfrag,1)
3406 allocate(nfrag(nlevel))
3408 ! Read the data pertaining to elementary fragments (level 1)
3409 call readi(controlcard,'NFRAG',nfrag(1),0)
3410 write(iout,*)"nfrag(1)",nfrag(1)
3411 call alloc_compar_arrays(nfrag(1),nlevel)
3413 call card_concat(controlcard,.true.)
3414 write (iout,*) controlcard(:ilen(controlcard))
3415 call readi(controlcard,'NPIECE',npiece(j,1),0)
3416 call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
3417 call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
3418 call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
3419 call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
3420 call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
3421 call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
3422 call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
3423 call readi(controlcard,'RMS',irms(j,1),0)
3424 call readi(controlcard,'LOCAL',iloc(j),1)
3425 call readi(controlcard,'ELCONT',ielecont(j,1),1)
3426 if (ielecont(j,1).eq.0) then
3427 call readi(controlcard,'SCCONT',isccont(j,1),1)
3429 ang_cut(j)=ang_cut(j)*deg2rad
3430 ang_cut1(j)=ang_cut1(j)*deg2rad
3432 call card_concat(controlcard,.true.)
3433 call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
3434 call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
3436 write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
3437 (ifrag(1,k,j),ifrag(2,k,j),&
3438 k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
3439 " ang_cut1",ang_cut1(j)*rad2deg
3440 write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
3441 write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
3442 write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
3443 " ilocal",iloc(j)," isccont",isccont(j,1)
3445 ! Read data pertaning to higher levels
3447 call card_concat(controlcard,.true.)
3448 call readi(controlcard,'NFRAG',NFRAG(i),0)
3449 write (iout,*) "i",i," nfrag",nfrag(i)
3451 call card_concat(controlcard,.true.)
3453 call readi(controlcard,'ELCONT',ielecont(j,i),0)
3454 if (ielecont(j,i).eq.0) then
3455 call readi(controlcard,'SCCONT',isccont(j,i),1)
3457 call readi(controlcard,'RMS',irms(j,i),0)
3463 call readi(controlcard,'NPIECE',npiece(j,i),0)
3464 call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
3465 call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
3466 call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
3468 call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
3469 call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
3470 write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
3471 n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
3472 " isccont",isccont(j,i)," irms",irms(j,i)
3473 write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
3474 write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
3475 write(iout,*)"nc_frac",nc_fragm(j,i),&
3476 " nc_req",nc_req_setf(j,i)
3479 if (binary) write (iout,*) "Classes written in binary format."
3482 call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
3483 call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
3484 call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
3485 call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
3486 call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
3487 call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
3488 call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
3489 call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
3490 call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
3491 call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
3492 call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
3493 call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
3494 call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
3495 call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
3496 call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
3497 call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
3498 call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
3499 call readi(controlcard,'RMS_SINGLE',irms_single,0)
3500 call readi(controlcard,'CONT_SINGLE',icont_single,1)
3501 call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
3502 call readi(controlcard,'RMS_PAIR',irms_pair,0)
3503 call readi(controlcard,'CONT_PAIR',icont_pair,1)
3504 call readi(controlcard,'SPLIT_BET',isplit_bet,0)
3505 angcut_hel=angcut_hel*deg2rad
3506 angcut1_hel=angcut1_hel*deg2rad
3507 angcut_bet=angcut_bet*deg2rad
3508 angcut1_bet=angcut1_bet*deg2rad
3509 angcut_strand=angcut_strand*deg2rad
3510 angcut1_strand=angcut1_strand*deg2rad
3511 write (iout,*) "Automatic detection of structural elements"
3512 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
3513 ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
3514 ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
3515 ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
3516 ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
3517 ' SPLIT_BET',isplit_bet
3518 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
3519 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
3520 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
3521 ' MAXANG_HEL',angcut1_hel*rad2deg
3522 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
3523 ' MAXANG_BET',angcut1_bet*rad2deg
3524 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
3525 ' MAXANG_STRAND',angcut1_strand*rad2deg
3526 write (iout,*) 'FRAC_MIN',frac_min_set
3528 end subroutine read_compar
3529 !--------------------------------------------------------------------------------
3531 !--------------------------------------------------------------------------------
3532 subroutine read_ref_structure(*)
3534 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
3537 use control_data, only:pdbref
3538 use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
3539 nstart_sup,nstart_seq,nperm,nres0
3540 use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
3541 use compare, only:seq_comp !,contact,elecont
3542 use geometry, only:chainbuild,dist
3543 use io_config, only:readpdb
3545 use conform_compar, only:contact,elecont
3547 ! include 'DIMENSIONS'
3548 ! include 'DIMENSIONS.ZSCOPT'
3549 ! include 'DIMENSIONS.COMPAR'
3550 ! include 'COMMON.IOUNITS'
3551 ! include 'COMMON.GEO'
3552 ! include 'COMMON.VAR'
3553 ! include 'COMMON.INTERACT'
3554 ! include 'COMMON.LOCAL'
3555 ! include 'COMMON.NAMES'
3556 ! include 'COMMON.CHAIN'
3557 ! include 'COMMON.FFIELD'
3558 ! include 'COMMON.SBRIDGE'
3559 ! include 'COMMON.HEADER'
3560 ! include 'COMMON.CONTROL'
3561 ! include 'COMMON.CONTACTS1'
3562 ! include 'COMMON.PEPTCONT'
3563 ! include 'COMMON.TIME1'
3564 ! include 'COMMON.COMPAR'
3565 character(len=4) :: sequence(nres)
3567 !el real(kind=8) :: x(maxvar)
3568 integer :: itype_pdb(nres,5)
3569 !el logical seq_comp
3570 integer :: i,j,k,nres_pdb,iaux,mnum
3571 real(kind=8) :: ddsc !el,dist
3572 integer :: kkk !,ilen
3576 write (iout,*) "pdbref",pdbref
3578 read(inp,'(a)') pdbfile
3579 write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
3580 pdbfile(:ilen(pdbfile))
3581 open(ipdbin,file=pdbfile,status='old',err=33)
3583 33 write (iout,'(a)') 'Error opening PDB file.'
3588 itype_pdb(i,mnum)=itype(i,mnum)
3594 iaux=itype_pdb(i,mnum)
3595 itype_pdb(i,mnum)=itype(i,mnum)
3603 if (nsup.le.(nct-nnt+1)) then
3604 do i=0,nct-nnt+1-nsup
3605 if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
3607 do j=nnt+nsup-1,nnt,-1
3609 cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
3612 do j=nnt+nsup-1,nnt,-1
3614 cref(k,j+i,kkk)=cref(k,j,kkk)
3616 write(iout,*) "J",j,"J+I",j+i
3617 phi_ref(j+i)=phi_ref(j)
3618 theta_ref(j+i)=theta_ref(j)
3619 alph_ref(j+i)=alph_ref(j)
3620 omeg_ref(j+i)=omeg_ref(j)
3624 write (iout,'(i5,3f10.5,5x,3f10.5)') &
3625 j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
3633 write (iout,'(a)') &
3634 'Error - sequences to be superposed do not match.'
3637 do i=0,nsup-(nct-nnt+1)
3638 if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
3641 nstart_sup=nstart_sup+i
3646 write (iout,'(a)') &
3647 'Error - sequences to be superposed do not match.'
3651 write (iout,'(a,i5)') &
3652 'Experimental structure begins at residue',nstart_seq
3654 call read_angles(inp,*38)
3656 38 write (iout,'(a)') 'Error reading reference structure.'
3665 cref(j,i,kkk)=c(j,i)
3669 nend_sup=nstart_sup+nsup-1
3672 c(j,i)=cref(j,i,kkk)
3678 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
3680 if (itype(i,mnum).ne.10) then
3681 ddsc = dist(i,nres+i)
3683 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
3687 dc_norm(j,nres+i)=0.0d0
3690 ! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
3691 ! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
3692 ! dc_norm(3,nres+i)**2
3694 dc(j,i)=c(j,i+1)-c(j,i)
3698 dc_norm(j,i)=dc(j,i)/ddsc
3701 ! print *,"Calling contact"
3702 call contact(.true.,ncont_ref,icont_ref(1,1),&
3703 nstart_sup,nend_sup)
3704 ! print *,"Calling elecont"
3705 call elecont(.true.,ncont_pept_ref,&
3706 icont_pept_ref(1,1),&
3707 nstart_sup,nend_sup)
3708 write (iout,'(a,i3,a,i3,a,i3,a)') &
3709 'Number of residues to be superposed:',nsup,&
3710 ' (from residue',nstart_sup,' to residue',&
3713 end subroutine read_ref_structure
3714 !--------------------------------------------------------------------------------
3716 !--------------------------------------------------------------------------------
3717 subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
3719 use geometry_data, only:nres,c
3720 use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
3721 ! implicit real*8 (a-h,o-z)
3722 ! include 'DIMENSIONS'
3723 ! include 'DIMENSIONS.ZSCOPT'
3724 ! include 'COMMON.CHAIN'
3725 ! include 'COMMON.INTERACT'
3726 ! include 'COMMON.NAMES'
3727 ! include 'COMMON.IOUNITS'
3728 ! include 'COMMON.HEADER'
3729 ! include 'COMMON.SBRIDGE'
3730 character(len=50) :: tytul
3731 character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',&
3732 'D','E','F','G','H','I','J','K','L','M','N','O',&
3733 'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
3734 integer,dimension(nres) :: ica !(maxres)
3735 real(kind=8) :: temp,efree,etot,entropy,rmsdev
3736 integer :: ii,i,j,iti,ires,iatom,ichain,mnum
3737 write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
3739 write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
3741 write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
3749 if (iti.eq.ntyp1) then
3750 if (itype(i-1,molnum(i-1)).eq.ntyp1) then
3753 write (ipdb,'(a)') 'TER'
3759 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
3763 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
3764 ires,(c(j,nres+i),j=1,3)
3768 write (ipdb,'(a)') 'TER'
3771 if (itype(i,mnum).eq.ntyp1) cycle
3772 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3773 write (ipdb,30) ica(i),ica(i+1)
3774 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3775 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
3776 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
3777 write (ipdb,30) ica(i),ica(i)+1
3780 if (itype(nct,molnum(nct)).ne.10) then
3781 write (ipdb,30) ica(nct),ica(nct)+1
3784 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
3786 write (ipdb,'(a6)') 'ENDMDL'
3787 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3788 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3789 30 FORMAT ('CONECT',8I5)
3791 end subroutine pdboutW
3793 !------------------------------------------------------------------------------
3795 !-----------------------------------------------------------------------------
3796 !-----------------------------------------------------------------------------