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
486 ! print *,"KURWA",ww(48)
487 ! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO "
488 ! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",&
489 ! "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",&
490 ! "WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",&
491 ! "WCATPROT ","WCATCAT
492 ! +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
493 ! +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
494 ! +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
495 ! +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
500 allocate(weights(n_ene))
515 weights(15)=0 !wstrain !
518 weights(18)=0 !scal14 !
522 weights(26)= wvdwpp_nucl
529 weights(32) =wbond_nucl
530 weights(33) =wang_nucl
532 weights(35) =wtor_nucl
533 weights(36) =wtor_d_nucl
534 weights(37) =wcorr_nucl
535 weights(38) =wcorr3_nucl
537 weights(42) =wcatprot
542 call card_concat(controlcard,.false.)
544 ! Return if not own parameters
546 if (iparm.ne.myparm .and. separate_parset) return
548 call reads(controlcard,"BONDPAR",bondname_t,bondname)
549 open (ibond,file=bondname_t,status='old')
551 call reads(controlcard,"THETPAR",thetname_t,thetname)
552 open (ithep,file=thetname_t,status='old')
554 call reads(controlcard,"ROTPAR",rotname_t,rotname)
555 open (irotam,file=rotname_t,status='old')
557 call reads(controlcard,"TORPAR",torname_t,torname)
558 open (itorp,file=torname_t,status='old')
560 call reads(controlcard,"TORDPAR",tordname_t,tordname)
561 open (itordp,file=tordname_t,status='old')
563 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
564 open (isccor,file=sccorname_t,status='old')
566 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
567 open (ifourier,file=fouriername_t,status='old')
569 call reads(controlcard,"ELEPAR",elename_t,elename)
570 open (ielep,file=elename_t,status='old')
572 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
573 open (isidep,file=sidename_t,status='old')
575 call reads(controlcard,"SCPPAR",scpname_t,scpname)
576 open (iscpp,file=scpname_t,status='old')
578 call getenv_loc('IONPAR',ionname)
579 open (iion,file=ionname,status='old')
581 write (iout,*) "Parameter set:",iparm
582 write (iout,*) "Energy-term weights:"
584 write (iout,'(a16,f10.5)') wname(i),ww(i)
586 write (iout,*) "Sidechain potential file : ",&
587 sidename_t(:ilen(sidename_t))
589 write (iout,*) "SCp potential file : ",&
590 scpname_t(:ilen(scpname_t))
592 write (iout,*) "Electrostatic potential file : ",&
593 elename_t(:ilen(elename_t))
594 write (iout,*) "Cumulant coefficient file : ",&
595 fouriername_t(:ilen(fouriername_t))
596 write (iout,*) "Torsional parameter file : ",&
597 torname_t(:ilen(torname_t))
598 write (iout,*) "Double torsional parameter file : ",&
599 tordname_t(:ilen(tordname_t))
600 write (iout,*) "Backbone-rotamer parameter file : ",&
601 sccorname(:ilen(sccorname))
602 write (iout,*) "Bond & inertia constant file : ",&
603 bondname_t(:ilen(bondname_t))
604 write (iout,*) "Bending parameter file : ",&
605 thetname_t(:ilen(thetname_t))
606 write (iout,*) "Rotamer parameter file : ",&
607 rotname_t(:ilen(rotname_t))
610 ! Read the virtual-bond parameters, masses, and moments of inertia
611 ! and Stokes' radii of the peptide group and side chains
613 allocate(dsc(ntyp1)) !(ntyp1)
614 allocate(dsc_inv(ntyp1)) !(ntyp1)
615 allocate(nbondterm(ntyp)) !(ntyp)
616 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
617 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
618 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
619 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
620 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
622 !el allocate(msc(ntyp+1)) !(ntyp+1)
623 !el allocate(isc(ntyp+1)) !(ntyp+1)
624 !el allocate(restok(ntyp+1)) !(ntyp+1)
625 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
628 read (ibond,*) vbldp0,akp
631 read (ibond,*) vbldsc0(1,i),aksc(1,i)
632 dsc(i) = vbldsc0(1,i)
636 dsc_inv(i)=1.0D0/dsc(i)
640 read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
642 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
644 dsc(i) = vbldsc0(1,i)
648 dsc_inv(i)=1.0D0/dsc(i)
653 write(iout,'(/a/)')"Force constants virtual bonds:"
654 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
656 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
658 write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
659 vbldsc0(1,i),aksc(1,i),abond0(1,i)
661 write (iout,'(13x,3f10.5)') &
662 vbldsc0(j,i),aksc(j,i),abond0(j,i)
666 if (.not. allocated(msc)) allocate(msc(ntyp1,5))
667 if (.not. allocated(restok)) allocate(restok(ntyp1,5))
670 read(iion,*) msc(i,5),restok(i,5)
671 print *,msc(i,5),restok(i,5)
675 read (iion,*) ncatprotparm
676 allocate(catprm(ncatprotparm,4))
678 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
681 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
684 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2)
685 ! dsc(i) = vbldsc0_nucl(1,i)
689 ! dsc_inv(i)=1.0D0/dsc(i)
694 !----------------------------------------------------
695 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
696 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
697 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
698 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
699 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
700 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
706 athet(j,i,ichir1,ichir2)=0.0D0
707 bthet(j,i,ichir1,ichir2)=0.0D0
721 !elwrite(iout,*) "parmread kontrol"
725 ! Read the parameters of the probability distribution/energy expression
726 ! of the virtual-bond valence angles theta
729 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
730 (bthet(j,i,1,1),j=1,2)
731 read (ithep,*) (polthet(j,i),j=0,3)
732 !elwrite(iout,*) "parmread kontrol in cryst_theta"
733 read (ithep,*) (gthet(j,i),j=1,3)
734 !elwrite(iout,*) "parmread kontrol in cryst_theta"
735 read (ithep,*) theta0(i),sig0(i),sigc0(i)
737 !elwrite(iout,*) "parmread kontrol in cryst_theta"
739 !elwrite(iout,*) "parmread kontrol in cryst_theta"
741 athet(1,i,1,-1)=athet(1,i,1,1)
742 athet(2,i,1,-1)=athet(2,i,1,1)
743 bthet(1,i,1,-1)=-bthet(1,i,1,1)
744 bthet(2,i,1,-1)=-bthet(2,i,1,1)
745 athet(1,i,-1,1)=-athet(1,i,1,1)
746 athet(2,i,-1,1)=-athet(2,i,1,1)
747 bthet(1,i,-1,1)=bthet(1,i,1,1)
748 bthet(2,i,-1,1)=bthet(2,i,1,1)
750 !elwrite(iout,*) "parmread kontrol in cryst_theta"
753 athet(1,i,-1,-1)=athet(1,-i,1,1)
754 athet(2,i,-1,-1)=-athet(2,-i,1,1)
755 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
756 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
757 athet(1,i,-1,1)=athet(1,-i,1,1)
758 athet(2,i,-1,1)=-athet(2,-i,1,1)
759 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
760 bthet(2,i,-1,1)=bthet(2,-i,1,1)
761 athet(1,i,1,-1)=-athet(1,-i,1,1)
762 athet(2,i,1,-1)=athet(2,-i,1,1)
763 bthet(1,i,1,-1)=bthet(1,-i,1,1)
764 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
769 polthet(j,i)=polthet(j,-i)
772 gthet(j,i)=gthet(j,-i)
775 !elwrite(iout,*) "parmread kontrol in cryst_theta"
777 !elwrite(iout,*) "parmread kontrol in cryst_theta"
780 ! & 'Parameters of the virtual-bond valence angles:'
781 ! write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
782 ! & ' ATHETA0 ',' A1 ',' A2 ',
785 ! write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
786 ! & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
788 ! write (iout,'(/a/9x,5a/79(1h-))')
789 ! & 'Parameters of the expression for sigma(theta_c):',
790 ! & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
791 ! & ' ALPH3 ',' SIGMA0C '
793 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
794 ! & (polthet(j,i),j=0,3),sigc0(i)
796 ! write (iout,'(/a/9x,5a/79(1h-))')
797 ! & 'Parameters of the second gaussian:',
798 ! & ' THETA0 ',' SIGMA0 ',' G1 ',
801 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
802 ! & sig0(i),(gthet(j,i),j=1,3)
805 'Parameters of the virtual-bond valence angles:'
806 write (iout,'(/a/9x,5a/79(1h-))') &
807 'Coefficients of expansion',&
808 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
809 ' b1*10^1 ',' b2*10^1 '
811 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
812 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
813 (10*bthet(j,i,1,1),j=1,2)
815 write (iout,'(/a/9x,5a/79(1h-))') &
816 'Parameters of the expression for sigma(theta_c):',&
817 ' alpha0 ',' alph1 ',' alph2 ',&
818 ' alhp3 ',' sigma0c '
820 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
821 (polthet(j,i),j=0,3),sigc0(i)
823 write (iout,'(/a/9x,5a/79(1h-))') &
824 'Parameters of the second gaussian:',&
825 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
828 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
829 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
834 ! Read the parameters of Utheta determined from ab initio surfaces
835 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
837 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
838 ! write (iout,*) "tu dochodze"
839 IF (tor_mode.eq.0) THEN
840 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
841 read (ithep,*) nthetyp,ntheterm,ntheterm2,&
842 ntheterm3,nsingle,ndouble
843 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
845 !----------------------------------------------------
846 ! allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
847 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
848 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
849 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
850 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
851 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
852 !(maxtheterm,-maxthetyp1:maxthetyp1,&
853 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
854 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
855 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
856 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
857 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
858 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
859 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
860 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
861 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
862 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
863 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
864 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
865 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
866 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
867 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
868 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
869 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
872 read (ithep,*) (ithetyp(i),i=1,ntyp1)
874 ithetyp(i)=-ithetyp(-i)
876 ! write (iout,*) "tu dochodze"
877 aa0thet(:,:,:,:)=0.0d0
878 aathet(:,:,:,:,:)=0.0d0
879 bbthet(:,:,:,:,:,:)=0.0d0
880 ccthet(:,:,:,:,:,:)=0.0d0
881 ddthet(:,:,:,:,:,:)=0.0d0
882 eethet(:,:,:,:,:,:)=0.0d0
883 ffthet(:,:,:,:,:,:,:)=0.0d0
884 ggthet(:,:,:,:,:,:,:)=0.0d0
888 do j=-nthetyp,nthetyp
889 do k=-nthetyp,nthetyp
890 read (ithep,'(6a)') res1
891 read (ithep,*) aa0thet(i,j,k,iblock)
892 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
894 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
895 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
896 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
897 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
900 (((ffthet(llll,lll,ll,i,j,k,iblock),&
901 ffthet(lll,llll,ll,i,j,k,iblock),&
902 ggthet(llll,lll,ll,i,j,k,iblock),&
903 ggthet(lll,llll,ll,i,j,k,iblock),&
904 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
909 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
910 ! coefficients of theta-and-gamma-dependent terms are zero.
915 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
916 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
918 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
919 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
922 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
924 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
927 ! Substitution for D aminoacids from symmetry.
930 do j=-nthetyp,nthetyp
931 do k=-nthetyp,nthetyp
932 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
934 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
938 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
939 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
940 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
941 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
947 ffthet(llll,lll,ll,i,j,k,iblock)= &
948 ffthet(llll,lll,ll,-i,-j,-k,iblock)
949 ffthet(lll,llll,ll,i,j,k,iblock)= &
950 ffthet(lll,llll,ll,-i,-j,-k,iblock)
951 ggthet(llll,lll,ll,i,j,k,iblock)= &
952 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
953 ggthet(lll,llll,ll,i,j,k,iblock)= &
954 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
964 ! Control printout of the coefficients of virtual-bond-angle potentials
968 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
972 write (iout,'(//4a)') &
973 'Type ',onelett(i),onelett(j),onelett(k)
974 write (iout,'(//a,10x,a)') " l","a[l]"
975 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
976 write (iout,'(i2,1pe15.5)') &
977 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
979 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
980 "b",l,"c",l,"d",l,"e",l
982 write (iout,'(i2,4(1pe15.5))') m,&
983 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
984 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
988 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
989 "f+",l,"f-",l,"g+",l,"g-",l
992 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
993 ffthet(n,m,l,i,j,k,iblock),&
994 ffthet(m,n,l,i,j,k,iblock),&
995 ggthet(n,m,l,i,j,k,iblock),&
996 ggthet(m,n,l,i,j,k,iblock)
1007 !C here will be the apropriate recalibrating for D-aminoacid
1008 read (ithep,*) nthetyp
1009 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1010 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1011 do i=-nthetyp+1,nthetyp-1
1012 read (ithep,*) nbend_kcc_Tb(i)
1013 do j=0,nbend_kcc_Tb(i)
1014 read (ithep,*) ijunk,v1bend_chyb(j,i)
1018 write (iout,'(a)') &
1019 "Parameters of the valence-only potentials"
1020 do i=-nthetyp+1,nthetyp-1
1021 write (iout,'(2a)') "Type ",toronelet(i)
1022 do k=0,nbend_kcc_Tb(i)
1023 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1030 !--------------- Reading theta parameters for nucleic acid-------
1031 read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,&
1032 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1033 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1034 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1035 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1036 nthetyp_nucl+1,nthetyp_nucl+1))
1037 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1038 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1039 nthetyp_nucl+1,nthetyp_nucl+1))
1040 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1041 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1042 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1043 nthetyp_nucl+1,nthetyp_nucl+1))
1044 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1045 nthetyp_nucl+1,nthetyp_nucl+1))
1046 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1047 nthetyp_nucl+1,nthetyp_nucl+1))
1048 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1049 nthetyp_nucl+1,nthetyp_nucl+1))
1050 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1051 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1052 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1053 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1054 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1055 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1057 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1058 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1060 read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1062 aa0thet_nucl(:,:,:)=0.0d0
1063 aathet_nucl(:,:,:,:)=0.0d0
1064 bbthet_nucl(:,:,:,:,:)=0.0d0
1065 ccthet_nucl(:,:,:,:,:)=0.0d0
1066 ddthet_nucl(:,:,:,:,:)=0.0d0
1067 eethet_nucl(:,:,:,:,:)=0.0d0
1068 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1069 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1074 read (ithep_nucl,'(3a)') t1,t2,t3
1075 read (ithep_nucl,*) aa0thet_nucl(i,j,k)
1076 read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1077 read (ithep_nucl,*) &
1078 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1079 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1080 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1081 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1082 read (ithep_nucl,*) &
1083 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1084 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1085 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1091 !-------------------------------------------
1092 allocate(nlob(ntyp1)) !(ntyp1)
1093 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1094 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1095 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1113 gaussc(l,k,j,i)=0.0D0
1121 ! Read the parameters of the probability distribution/energy expression
1122 ! of the side chains.
1125 !c write (iout,*) "tu dochodze",i
1126 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
1130 dsc_inv(i)=1.0D0/dsc(i)
1141 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
1142 censc(1,1,-i)=censc(1,1,i)
1143 censc(2,1,-i)=censc(2,1,i)
1144 censc(3,1,-i)=-censc(3,1,i)
1146 read (irotam,*) bsc(j,i)
1147 read (irotam,*) (censc(k,j,i),k=1,3),&
1148 ((blower(k,l,j),l=1,k),k=1,3)
1149 censc(1,j,-i)=censc(1,j,i)
1150 censc(2,j,-i)=censc(2,j,i)
1151 censc(3,j,-i)=-censc(3,j,i)
1152 ! BSC is amplitude of Gaussian
1159 akl=akl+blower(k,m,j)*blower(l,m,j)
1163 if (((k.eq.3).and.(l.ne.3)) &
1164 .or.((l.eq.3).and.(k.ne.3))) then
1165 gaussc(k,l,j,-i)=-akl
1166 gaussc(l,k,j,-i)=-akl
1168 gaussc(k,l,j,-i)=akl
1169 gaussc(l,k,j,-i)=akl
1178 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1181 if (nlobi.gt.0) then
1182 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1183 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1184 ! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1185 ! write (iout,'(a,f10.4,4(16x,f10.4))')
1186 ! & 'Center ',(bsc(j,i),j=1,nlobi)
1187 ! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
1188 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1189 'log h',(bsc(j,i),j=1,nlobi)
1190 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1191 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1192 ! write (iout,'(a)')
1198 ! blower(k,l,j)=gaussc(ind,j,i)
1203 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1204 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1211 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1212 ! added by Urszula Kozlowska 07/11/2007
1214 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1222 read(irotam,*) sc_parmin(j,i)
1228 !---------reading nucleic acid parameters for rotamers-------------------
1229 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1230 do i=1,ntyp_molec(2)
1231 read (irotam_nucl,*)
1233 read(irotam_nucl,*) sc_parmin_nucl(j,i)
1239 write (iout,*) "Base rotamer parameters"
1240 do i=1,ntyp_molec(2)
1241 write (iout,'(a)') restyp(i,2)
1242 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1247 read (ifourier,*) nloctyp
1248 !el write(iout,*)"nloctyp",nloctyp
1249 SPLIT_FOURIERTOR = nloctyp.lt.0
1250 nloctyp = iabs(nloctyp)
1252 if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1))
1253 print *,"shape",shape(itype2loc)
1254 allocate(iloctyp(-nloctyp:nloctyp))
1255 allocate(bnew1(3,2,-nloctyp:nloctyp))
1256 allocate(bnew2(3,2,-nloctyp:nloctyp))
1257 allocate(ccnew(3,2,-nloctyp:nloctyp))
1258 allocate(ddnew(3,2,-nloctyp:nloctyp))
1259 allocate(e0new(3,-nloctyp:nloctyp))
1260 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1261 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1262 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1263 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1264 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1265 allocate(e0newtor(3,-nloctyp:nloctyp))
1266 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1268 read (ifourier,*) (itype2loc(i),i=1,ntyp)
1269 read (ifourier,*) (iloctyp(i),i=0,nloctyp-1)
1270 itype2loc(ntyp1)=nloctyp
1271 iloctyp(nloctyp)=ntyp1
1273 itype2loc(-i)=-itype2loc(i)
1276 allocate(iloctyp(-nloctyp:nloctyp))
1277 allocate(itype2loc(-ntyp1:ntyp1))
1284 iloctyp(-i)=-iloctyp(i)
1286 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1287 !c write (iout,*) "nloctyp",nloctyp,
1288 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1289 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1290 !c write (iout,*) "nloctyp",nloctyp,
1291 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1294 !c write (iout,*) "NEWCORR",i
1298 read (ifourier,*) bnew1(ii,j,i)
1301 !c write (iout,*) "NEWCORR BNEW1"
1302 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1305 read (ifourier,*) bnew2(ii,j,i)
1308 !c write (iout,*) "NEWCORR BNEW2"
1309 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1311 read (ifourier,*) ccnew(kk,1,i)
1312 read (ifourier,*) ccnew(kk,2,i)
1314 !c write (iout,*) "NEWCORR CCNEW"
1315 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1317 read (ifourier,*) ddnew(kk,1,i)
1318 read (ifourier,*) ddnew(kk,2,i)
1320 !c write (iout,*) "NEWCORR DDNEW"
1321 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1325 read (ifourier,*) eenew(ii,jj,kk,i)
1329 !c write (iout,*) "NEWCORR EENEW1"
1330 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1332 read (ifourier,*) e0new(ii,i)
1334 !c write (iout,*) (e0new(ii,i),ii=1,3)
1336 !c write (iout,*) "NEWCORR EENEW"
1339 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1340 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1341 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1342 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1347 bnew1(ii,1,-i)= bnew1(ii,1,i)
1348 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1349 bnew2(ii,1,-i)= bnew2(ii,1,i)
1350 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1353 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1354 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1355 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1356 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1357 ccnew(ii,1,-i)=ccnew(ii,1,i)
1358 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1359 ddnew(ii,1,-i)=ddnew(ii,1,i)
1360 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1362 e0new(1,-i)= e0new(1,i)
1363 e0new(2,-i)=-e0new(2,i)
1364 e0new(3,-i)=-e0new(3,i)
1366 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1367 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1368 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1369 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1373 write (iout,'(a)') "Coefficients of the multibody terms"
1374 do i=-nloctyp+1,nloctyp-1
1375 write (iout,*) "Type: ",onelet(iloctyp(i))
1376 write (iout,*) "Coefficients of the expansion of B1"
1378 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1380 write (iout,*) "Coefficients of the expansion of B2"
1382 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1384 write (iout,*) "Coefficients of the expansion of C"
1385 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1386 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1387 write (iout,*) "Coefficients of the expansion of D"
1388 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1389 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1390 write (iout,*) "Coefficients of the expansion of E"
1391 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1394 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1399 IF (SPLIT_FOURIERTOR) THEN
1401 !c write (iout,*) "NEWCORR TOR",i
1405 read (ifourier,*) bnew1tor(ii,j,i)
1408 !c write (iout,*) "NEWCORR BNEW1 TOR"
1409 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1412 read (ifourier,*) bnew2tor(ii,j,i)
1415 !c write (iout,*) "NEWCORR BNEW2 TOR"
1416 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1418 read (ifourier,*) ccnewtor(kk,1,i)
1419 read (ifourier,*) ccnewtor(kk,2,i)
1421 !c write (iout,*) "NEWCORR CCNEW TOR"
1422 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1424 read (ifourier,*) ddnewtor(kk,1,i)
1425 read (ifourier,*) ddnewtor(kk,2,i)
1427 !c write (iout,*) "NEWCORR DDNEW TOR"
1428 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1432 read (ifourier,*) eenewtor(ii,jj,kk,i)
1436 !c write (iout,*) "NEWCORR EENEW1 TOR"
1437 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1439 read (ifourier,*) e0newtor(ii,i)
1441 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1443 !c write (iout,*) "NEWCORR EENEW TOR"
1446 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1447 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1448 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1449 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1454 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1455 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1456 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1457 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1460 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1461 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1462 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1463 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1465 e0newtor(1,-i)= e0newtor(1,i)
1466 e0newtor(2,-i)=-e0newtor(2,i)
1467 e0newtor(3,-i)=-e0newtor(3,i)
1469 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1470 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1471 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1472 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1476 write (iout,'(a)') &
1477 "Single-body coefficients of the torsional potentials"
1478 do i=-nloctyp+1,nloctyp-1
1479 write (iout,*) "Type: ",onelet(iloctyp(i))
1480 write (iout,*) "Coefficients of the expansion of B1tor"
1482 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1483 j,(bnew1tor(k,j,i),k=1,3)
1485 write (iout,*) "Coefficients of the expansion of B2tor"
1487 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1488 j,(bnew2tor(k,j,i),k=1,3)
1490 write (iout,*) "Coefficients of the expansion of Ctor"
1491 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1492 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1493 write (iout,*) "Coefficients of the expansion of Dtor"
1494 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1495 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1496 write (iout,*) "Coefficients of the expansion of Etor"
1497 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1500 write (iout,'(1hE,2i1,2f10.5)') &
1501 j,k,(eenewtor(l,j,k,i),l=1,2)
1507 do i=-nloctyp+1,nloctyp-1
1510 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1515 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1519 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1520 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1521 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1522 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1525 ENDIF !SPLIT_FOURIER_TOR
1527 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1528 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1529 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1530 allocate(b(13,-nloctyp-1:nloctyp+1))
1532 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1535 read (ifourier,*) (b(ii,i),ii=1,13)
1537 write (iout,*) 'Type ',onelet(iloctyp(i))
1538 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1546 CCold(1,1,i)= b(7,i)
1547 CCold(2,2,i)=-b(7,i)
1548 CCold(2,1,i)= b(9,i)
1549 CCold(1,2,i)= b(9,i)
1550 CCold(1,1,-i)= b(7,i)
1551 CCold(2,2,-i)=-b(7,i)
1552 CCold(2,1,-i)=-b(9,i)
1553 CCold(1,2,-i)=-b(9,i)
1554 DDold(1,1,i)= b(6,i)
1555 DDold(2,2,i)=-b(6,i)
1556 DDold(2,1,i)= b(8,i)
1557 DDold(1,2,i)= b(8,i)
1558 DDold(1,1,-i)= b(6,i)
1559 DDold(2,2,-i)=-b(6,i)
1560 DDold(2,1,-i)=-b(8,i)
1561 DDold(1,2,-i)=-b(8,i)
1562 EEold(1,1,i)= b(10,i)+b(11,i)
1563 EEold(2,2,i)=-b(10,i)+b(11,i)
1564 EEold(2,1,i)= b(12,i)-b(13,i)
1565 EEold(1,2,i)= b(12,i)+b(13,i)
1566 EEold(1,1,-i)= b(10,i)+b(11,i)
1567 EEold(2,2,-i)=-b(10,i)+b(11,i)
1568 EEold(2,1,-i)=-b(12,i)+b(13,i)
1569 EEold(1,2,-i)=-b(12,i)-b(13,i)
1570 write(iout,*) "TU DOCHODZE"
1576 "Coefficients of the cumulants (independent of valence angles)"
1577 do i=-nloctyp+1,nloctyp-1
1578 write (iout,*) 'Type ',onelet(iloctyp(i))
1580 write(iout,'(2f10.5)') B(3,i),B(5,i)
1582 write(iout,'(2f10.5)') B(2,i),B(4,i)
1585 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1589 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1593 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1600 ! Read torsional parameters in old format
1602 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1604 read (itorp,*) ntortyp,nterm_old
1605 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1606 read (itorp,*) (itortyp(i),i=1,ntyp)
1608 !el from energy module--------
1609 allocate(v1(nterm_old,ntortyp,ntortyp))
1610 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1611 !el---------------------------
1617 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1623 write (iout,'(/a/)') 'Torsional constants:'
1626 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1627 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1635 ! Read torsional parameters
1637 IF (TOR_MODE.eq.0) THEN
1639 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1641 read (itorp,*) ntortyp
1642 read (itorp,*) (itortyp(i),i=1,ntyp)
1643 write (iout,*) 'ntortyp',ntortyp
1645 !el from energy module---------
1646 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1647 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1649 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1650 allocate(vlor2(maxlor,ntortyp,ntortyp))
1651 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1652 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1654 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1655 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1656 !el---------------------------
1658 do i=-ntortyp,ntortyp
1659 do j=-ntortyp,ntortyp
1665 !el---------------------------
1669 itortyp(i)=-itortyp(-i)
1671 ! write (iout,*) 'ntortyp',ntortyp
1673 do j=-ntortyp+1,ntortyp-1
1674 read (itorp,*) nterm(i,j,iblock),&
1676 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1677 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1680 do k=1,nterm(i,j,iblock)
1681 read (itorp,*) kk,v1(k,i,j,iblock),&
1683 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1684 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1685 v0ij=v0ij+si*v1(k,i,j,iblock)
1688 do k=1,nlor(i,j,iblock)
1689 read (itorp,*) kk,vlor1(k,i,j),&
1690 vlor2(k,i,j),vlor3(k,i,j)
1691 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1694 v0(-i,-j,iblock)=v0ij
1701 write (iout,'(/a/)') 'Torsional constants:'
1704 write (iout,*) 'ityp',i,' jtyp',j
1705 write (iout,*) 'Fourier constants'
1706 do k=1,nterm(i,j,iblock)
1707 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1710 write (iout,*) 'Lorenz constants'
1711 do k=1,nlor(i,j,iblock)
1712 write (iout,'(3(1pe15.5))') &
1713 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1720 ! 6/23/01 Read parameters for double torsionals
1722 !el from energy module------------
1723 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1724 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1725 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1726 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1727 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1728 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1729 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1730 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1731 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1732 !---------------------------------
1736 do j=-ntortyp+1,ntortyp-1
1737 do k=-ntortyp+1,ntortyp-1
1738 read (itordp,'(3a1)') t1,t2,t3
1739 ! write (iout,*) "OK onelett",
1742 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1743 .or. t3.ne.toronelet(k)) then
1744 write (iout,*) "Error in double torsional parameter file",&
1747 call MPI_Finalize(Ierror)
1749 stop "Error in double torsional parameter file"
1751 read (itordp,*) ntermd_1(i,j,k,iblock),&
1752 ntermd_2(i,j,k,iblock)
1753 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1754 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1755 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1756 ntermd_1(i,j,k,iblock))
1757 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1758 ntermd_1(i,j,k,iblock))
1759 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1760 ntermd_1(i,j,k,iblock))
1761 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1762 ntermd_1(i,j,k,iblock))
1763 ! Martix of D parameters for one dimesional foureir series
1764 do l=1,ntermd_1(i,j,k,iblock)
1765 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1766 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1767 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1768 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1769 ! write(iout,*) "whcodze" ,
1770 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1772 read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1773 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1774 v2s(m,l,i,j,k,iblock),&
1775 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1776 ! Martix of D parameters for two dimesional fourier series
1777 do l=1,ntermd_2(i,j,k,iblock)
1779 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1780 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1781 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1782 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1791 write (iout,*) 'Constants for double torsionals'
1794 do j=-ntortyp+1,ntortyp-1
1795 do k=-ntortyp+1,ntortyp-1
1796 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1797 ' nsingle',ntermd_1(i,j,k,iblock),&
1798 ' ndouble',ntermd_2(i,j,k,iblock)
1800 write (iout,*) 'Single angles:'
1801 do l=1,ntermd_1(i,j,k,iblock)
1802 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1803 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1804 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1805 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1808 write (iout,*) 'Pairs of angles:'
1809 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1810 do l=1,ntermd_2(i,j,k,iblock)
1811 write (iout,'(i5,20f10.5)') &
1812 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1815 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1816 do l=1,ntermd_2(i,j,k,iblock)
1817 write (iout,'(i5,20f10.5)') &
1818 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1819 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1829 itype2loc(i)=itortyp(i)
1832 ELSE IF (TOR_MODE.eq.1) THEN
1834 !C read valence-torsional parameters
1835 read (itorp,*) ntortyp
1837 write (iout,*) "Valence-torsional parameters read in ntortyp",&
1839 read (itorp,*) (itortyp(i),i=1,ntyp)
1840 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1843 itype2loc(i)=itortyp(i)
1847 itortyp(i)=-itortyp(-i)
1849 do i=-ntortyp+1,ntortyp-1
1850 do j=-ntortyp+1,ntortyp-1
1851 !C first we read the cos and sin gamma parameters
1852 read (itorp,'(13x,a)') string
1853 write (iout,*) i,j,string
1855 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1856 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1857 do k=1,nterm_kcc(j,i)
1858 do l=1,nterm_kcc_Tb(j,i)
1859 do ll=1,nterm_kcc_Tb(j,i)
1860 read (itorp,*) ii,jj,kk, &
1861 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1869 !c AL 4/8/16: Calculate coefficients from one-body parameters
1871 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1872 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
1873 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
1874 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1875 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1878 print *,i,itortyp(i)
1879 itortyp(i)=itype2loc(i)
1882 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
1884 do i=-ntortyp+1,ntortyp-1
1885 do j=-ntortyp+1,ntortyp-1
1888 do k=1,nterm_kcc_Tb(j,i)
1889 do l=1,nterm_kcc_Tb(j,i)
1890 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
1891 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1892 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
1893 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1896 do k=1,nterm_kcc_Tb(j,i)
1897 do l=1,nterm_kcc_Tb(j,i)
1899 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1900 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1901 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1902 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1904 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1905 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1906 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1907 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1911 !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)
1915 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1919 if (tor_mode.gt.0 .and. lprint) then
1920 !c Print valence-torsional parameters
1921 write (iout,'(a)') &
1922 "Parameters of the valence-torsional potentials"
1923 do i=-ntortyp+1,ntortyp-1
1924 do j=-ntortyp+1,ntortyp-1
1925 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1926 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1927 do k=1,nterm_kcc(j,i)
1928 do l=1,nterm_kcc_Tb(j,i)
1929 do ll=1,nterm_kcc_Tb(j,i)
1930 write (iout,'(3i5,2f15.4)')&
1931 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1940 !elwrite(iout,*) "parmread kontrol sc-bb"
1941 ! Read of Side-chain backbone correlation parameters
1942 ! Modified 11 May 2012 by Adasko
1945 read (isccor,*) nsccortyp
1948 !c maxinter is maximum interaction sites
1949 !write(iout,*)"maxterm_sccor",maxterm_sccor
1950 !el from module energy-------------
1951 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1952 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1953 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1954 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1955 !-----------------------------------
1956 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1957 !-----------------------------------
1958 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1959 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1960 -nsccortyp:nsccortyp))
1961 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1962 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1963 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1964 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1965 !-----------------------------------
1966 do i=-nsccortyp,nsccortyp
1967 do j=-nsccortyp,nsccortyp
1971 !-----------------------------------
1973 read (isccor,*) (isccortyp(i),i=1,ntyp)
1975 isccortyp(i)=-isccortyp(-i)
1977 iscprol=isccortyp(20)
1978 ! write (iout,*) 'ntortyp',ntortyp
1980 !c maxinter is maximum interaction sites
1985 nterm_sccor(i,j),nlor_sccor(i,j)
1991 nterm_sccor(-i,j)=nterm_sccor(i,j)
1992 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1993 nterm_sccor(i,-j)=nterm_sccor(i,j)
1994 do k=1,nterm_sccor(i,j)
1995 read (isccor,*) kk,v1sccor(k,l,i,j),&
1997 if (j.eq.iscprol) then
1998 if (i.eq.isccortyp(10)) then
1999 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2000 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2002 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2003 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2004 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2005 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2006 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2007 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2008 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2009 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2012 if (i.eq.isccortyp(10)) then
2013 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2014 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2016 if (j.eq.isccortyp(10)) then
2017 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2018 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2020 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2021 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2022 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2023 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2024 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2025 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2029 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2030 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2031 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2032 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2035 do k=1,nlor_sccor(i,j)
2036 read (isccor,*) kk,vlor1sccor(k,i,j),&
2037 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2038 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2039 (1+vlor3sccor(k,i,j)**2)
2041 v0sccor(l,i,j)=v0ijsccor
2042 v0sccor(l,-i,j)=v0ijsccor1
2043 v0sccor(l,i,-j)=v0ijsccor2
2044 v0sccor(l,-i,-j)=v0ijsccor3
2050 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
2053 write (iout,*) 'ityp',i,' jtyp',j
2054 write (iout,*) 'Fourier constants'
2055 do k=1,nterm_sccor(i,j)
2056 write (iout,'(2(1pe15.5))') &
2057 (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
2059 write (iout,*) 'Lorenz constants'
2060 do k=1,nlor_sccor(i,j)
2061 write (iout,'(3(1pe15.5))') &
2062 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2069 ! Read electrostatic-interaction parameters
2072 write (iout,'(/a)') 'Electrostatic interaction constants:'
2073 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2074 'IT','JT','APP','BPP','AEL6','AEL3'
2076 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
2077 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
2078 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
2079 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
2084 app (i,j)=epp(i,j)*rri*rri
2085 bpp (i,j)=-2.0D0*epp(i,j)*rri
2086 ael6(i,j)=elpp6(i,j)*4.2D0**6
2087 ael3(i,j)=elpp3(i,j)*4.2D0**3
2088 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2093 ! Read side-chain interaction parameters.
2095 !el from module energy - COMMON.INTERACT-------
2096 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2097 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2098 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2099 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2100 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2101 allocate(epslip(ntyp,ntyp))
2112 !--------------------------------
2114 read (isidep,*) ipot,expon
2115 !el if (ipot.lt.1 .or. ipot.gt.5) then
2116 ! write (iout,'(2a)') 'Error while reading SC interaction',&
2117 ! 'potential file - unknown potential type.'
2121 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2122 ', exponents are ',expon,2*expon
2123 ! goto (10,20,30,30,40) ipot
2125 !----------------------- LJ potential ---------------------------------
2127 ! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2128 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2130 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2131 write (iout,'(a/)') 'The epsilon array:'
2132 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2133 write (iout,'(/a)') 'One-body parameters:'
2134 write (iout,'(a,4x,a)') 'residue','sigma'
2135 write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
2138 !----------------------- LJK potential --------------------------------
2140 ! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2141 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2142 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2144 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2145 write (iout,'(a/)') 'The epsilon array:'
2146 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2147 write (iout,'(/a)') 'One-body parameters:'
2148 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2149 write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2153 !---------------------- GB or BP potential -----------------------------
2156 if (scelemode.eq.0) then
2158 read (isidep,*)(eps(i,j),j=i,ntyp)
2160 read (isidep,*)(sigma0(i),i=1,ntyp)
2161 read (isidep,*)(sigii(i),i=1,ntyp)
2162 read (isidep,*)(chip(i),i=1,ntyp)
2163 read (isidep,*)(alp(i),i=1,ntyp)
2165 read (isidep,*)(epslip(i,j),j=i,ntyp)
2167 ! For the GB potential convert sigma'**2 into chi'
2170 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2174 write (iout,'(/a/)') 'Parameters of the BP potential:'
2175 write (iout,'(a/)') 'The epsilon array:'
2176 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2177 write (iout,'(/a)') 'One-body parameters:'
2178 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2180 write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
2181 chip(i),alp(i),i=1,ntyp)
2184 allocate(icharge(ntyp1))
2185 ! print *,ntyp,icharge(i)
2187 read (isidep,*) (icharge(i),i=1,ntyp)
2188 print *,ntyp,icharge(i)
2189 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2190 write (2,*) "icharge",(icharge(i),i=1,ntyp)
2191 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2192 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2193 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2194 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2195 allocate(epsintab(ntyp,ntyp))
2196 allocate(dtail(2,ntyp,ntyp))
2197 print *,"control line 1"
2198 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2199 allocate(wqdip(2,ntyp,ntyp))
2200 allocate(wstate(4,ntyp,ntyp))
2201 allocate(dhead(2,2,ntyp,ntyp))
2202 allocate(nstate(ntyp,ntyp))
2203 allocate(debaykap(ntyp,ntyp))
2204 print *,"control line 2"
2205 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2206 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2210 ! write (*,*) "Im in ALAB", i, " ", j
2212 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2213 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2214 chis(i,j),chis(j,i), &
2215 nstate(i,j),(wstate(k,i,j),k=1,4), &
2216 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2217 dtail(1,i,j),dtail(2,i,j), &
2218 epshead(i,j),sig0head(i,j), &
2219 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2220 alphapol(i,j),alphapol(j,i), &
2221 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
2222 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2228 sigma(i,j) = sigma(j,i)
2229 sigmap1(i,j)=sigmap1(j,i)
2230 sigmap2(i,j)=sigmap2(j,i)
2231 sigiso1(i,j)=sigiso1(j,i)
2232 sigiso2(i,j)=sigiso2(j,i)
2233 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2234 nstate(i,j) = nstate(j,i)
2235 dtail(1,i,j) = dtail(1,j,i)
2236 dtail(2,i,j) = dtail(2,j,i)
2238 alphasur(k,i,j) = alphasur(k,j,i)
2239 wstate(k,i,j) = wstate(k,j,i)
2240 alphiso(k,i,j) = alphiso(k,j,i)
2243 dhead(2,1,i,j) = dhead(1,1,j,i)
2244 dhead(2,2,i,j) = dhead(1,2,j,i)
2245 dhead(1,1,i,j) = dhead(2,1,j,i)
2246 dhead(1,2,i,j) = dhead(2,2,j,i)
2248 epshead(i,j) = epshead(j,i)
2249 sig0head(i,j) = sig0head(j,i)
2252 wqdip(k,i,j) = wqdip(k,j,i)
2255 wquad(i,j) = wquad(j,i)
2256 epsintab(i,j) = epsintab(j,i)
2257 debaykap(i,j)=debaykap(j,i)
2258 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2264 !--------------------- GBV potential -----------------------------------
2266 ! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2267 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2268 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2269 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2271 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2272 write (iout,'(a/)') 'The epsilon array:'
2273 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2274 write (iout,'(/a)') 'One-body parameters:'
2275 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2276 's||/s_|_^2',' chip ',' alph '
2277 write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2278 sigii(i),chip(i),alp(i),i=1,ntyp)
2281 write (iout,'(2a)') 'Error while reading SC interaction',&
2282 'potential file - unknown potential type.'
2288 !-----------------------------------------------------------------------
2289 ! Calculate the "working" parameters of SC interactions.
2291 !el from module energy - COMMON.INTERACT-------
2292 ! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2293 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1))
2294 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) !(ntyp,ntyp)
2295 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2296 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2297 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2298 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2300 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2308 if (scelemode.eq.0) then
2315 !--------------------------------
2320 epslip(i,j)=epslip(j,i)
2323 if (scelemode.eq.0) then
2326 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2327 sigma(j,i)=sigma(i,j)
2328 rs0(i,j)=dwa16*sigma(i,j)
2333 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2334 'Working parameters of the SC interactions:',&
2335 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2340 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2350 sigeps=dsign(1.0D0,epsij)
2352 aa_aq(i,j)=epsij*rrij*rrij
2353 bb_aq(i,j)=-sigeps*epsij*rrij
2354 aa_aq(j,i)=aa_aq(i,j)
2355 bb_aq(j,i)=bb_aq(i,j)
2356 epsijlip=epslip(i,j)
2357 sigeps=dsign(1.0D0,epsijlip)
2358 epsijlip=dabs(epsijlip)
2359 aa_lip(i,j)=epsijlip*rrij*rrij
2360 bb_lip(i,j)=-sigeps*epsijlip*rrij
2361 aa_lip(j,i)=aa_lip(i,j)
2362 bb_lip(j,i)=bb_lip(i,j)
2363 if ((ipot.gt.2).and. (scelemode.eq.0))then
2364 sigt1sq=sigma0(i)**2
2365 sigt2sq=sigma0(j)**2
2368 ratsig1=sigt2sq/sigt1sq
2369 ratsig2=1.0D0/ratsig1
2370 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2371 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2372 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2376 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2377 sigmaii(i,j)=rsum_max
2378 sigmaii(j,i)=rsum_max
2380 ! sigmaii(i,j)=r0(i,j)
2381 ! sigmaii(j,i)=r0(i,j)
2383 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2384 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2385 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2386 augm(i,j)=epsij*r_augm**(2*expon)
2387 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2394 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2395 restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2396 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2400 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2401 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2402 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2403 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2404 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2405 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2406 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2407 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2408 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2409 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2410 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2411 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2412 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2413 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2414 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2415 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2424 read (isidep_nucl,*) ipot_nucl
2425 ! print *,"TU?!",ipot_nucl
2426 if (ipot_nucl.eq.1) then
2427 do i=1,ntyp_molec(2)
2428 do j=i,ntyp_molec(2)
2429 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2430 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2434 do i=1,ntyp_molec(2)
2435 do j=i,ntyp_molec(2)
2436 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2437 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2438 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2442 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2443 do i=1,ntyp_molec(2)
2444 do j=i,ntyp_molec(2)
2445 rrij=sigma_nucl(i,j)
2449 epsij=4*eps_nucl(i,j)
2450 sigeps=dsign(1.0D0,epsij)
2452 aa_nucl(i,j)=epsij*rrij*rrij
2453 bb_nucl(i,j)=-sigeps*epsij*rrij
2454 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2455 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2456 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2457 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2458 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2459 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2462 aa_nucl(i,j)=aa_nucl(j,i)
2463 bb_nucl(i,j)=bb_nucl(j,i)
2464 ael3_nucl(i,j)=ael3_nucl(j,i)
2465 ael6_nucl(i,j)=ael6_nucl(j,i)
2466 ael63_nucl(i,j)=ael63_nucl(j,i)
2467 ael32_nucl(i,j)=ael32_nucl(j,i)
2468 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2469 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2470 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2471 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2472 eps_nucl(i,j)=eps_nucl(j,i)
2473 sigma_nucl(i,j)=sigma_nucl(j,i)
2474 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2478 write(iout,*) "tube param"
2479 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2480 ccavtubpep,dcavtubpep,tubetranenepep
2481 sigmapeptube=sigmapeptube**6
2482 sigeps=dsign(1.0D0,epspeptube)
2483 epspeptube=dabs(epspeptube)
2484 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2485 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2486 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2488 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2489 ccavtub(i),dcavtub(i),tubetranene(i)
2490 sigmasctube=sigmasctube**6
2491 sigeps=dsign(1.0D0,epssctube)
2492 epssctube=dabs(epssctube)
2493 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2494 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2495 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2497 !-----------------READING SC BASE POTENTIALS-----------------------------
2498 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2499 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2500 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2501 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2502 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2503 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2504 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2505 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2506 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2507 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2508 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2509 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2510 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2511 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2512 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2513 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2516 do i=1,ntyp_molec(1)
2517 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2518 write (*,*) "Im in ", i, " ", j
2519 read(isidep_scbase,*) &
2520 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2521 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2522 write(*,*) "eps",eps_scbase(i,j)
2523 read(isidep_scbase,*) &
2524 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2525 chis_scbase(i,j,1),chis_scbase(i,j,2)
2526 read(isidep_scbase,*) &
2527 dhead_scbasei(i,j), &
2528 dhead_scbasej(i,j), &
2529 rborn_scbasei(i,j),rborn_scbasej(i,j)
2530 read(isidep_scbase,*) &
2531 (wdipdip_scbase(k,i,j),k=1,3), &
2532 (wqdip_scbase(k,i,j),k=1,2)
2533 read(isidep_scbase,*) &
2534 alphapol_scbase(i,j), &
2535 epsintab_scbase(i,j)
2538 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2539 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2541 do i=1,ntyp_molec(1)
2542 do j=1,ntyp_molec(2)-1
2543 epsij=eps_scbase(i,j)
2544 rrij=sigma_scbase(i,j)
2549 sigeps=dsign(1.0D0,epsij)
2551 aa_scbase(i,j)=epsij*rrij*rrij
2552 bb_scbase(i,j)=-sigeps*epsij*rrij
2555 !-----------------READING PEP BASE POTENTIALS-------------------
2556 allocate(eps_pepbase(ntyp_molec(2)))
2557 allocate(sigma_pepbase(ntyp_molec(2)))
2558 allocate(chi_pepbase(ntyp_molec(2),2))
2559 allocate(chipp_pepbase(ntyp_molec(2),2))
2560 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2561 allocate(sigmap1_pepbase(ntyp_molec(2)))
2562 allocate(sigmap2_pepbase(ntyp_molec(2)))
2563 allocate(chis_pepbase(ntyp_molec(2),2))
2564 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2567 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2568 write (*,*) "Im in ", i, " ", j
2569 read(isidep_pepbase,*) &
2570 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2571 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2572 write(*,*) "eps",eps_pepbase(j)
2573 read(isidep_pepbase,*) &
2574 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2575 chis_pepbase(j,1),chis_pepbase(j,2)
2576 read(isidep_pepbase,*) &
2577 (wdipdip_pepbase(k,j),k=1,3)
2579 allocate(aa_pepbase(ntyp_molec(2)))
2580 allocate(bb_pepbase(ntyp_molec(2)))
2582 do j=1,ntyp_molec(2)-1
2583 epsij=eps_pepbase(j)
2584 rrij=sigma_pepbase(j)
2589 sigeps=dsign(1.0D0,epsij)
2591 aa_pepbase(j)=epsij*rrij*rrij
2592 bb_pepbase(j)=-sigeps*epsij*rrij
2594 !--------------READING SC PHOSPHATE-------------------------------------
2595 !--------------READING SC PHOSPHATE-------------------------------------
2596 allocate(eps_scpho(ntyp_molec(1)))
2597 allocate(sigma_scpho(ntyp_molec(1)))
2598 allocate(chi_scpho(ntyp_molec(1),2))
2599 allocate(chipp_scpho(ntyp_molec(1),2))
2600 allocate(alphasur_scpho(4,ntyp_molec(1)))
2601 allocate(sigmap1_scpho(ntyp_molec(1)))
2602 allocate(sigmap2_scpho(ntyp_molec(1)))
2603 allocate(chis_scpho(ntyp_molec(1),2))
2604 allocate(wqq_scpho(ntyp_molec(1)))
2605 allocate(wqdip_scpho(2,ntyp_molec(1)))
2606 allocate(alphapol_scpho(ntyp_molec(1)))
2607 allocate(epsintab_scpho(ntyp_molec(1)))
2608 allocate(dhead_scphoi(ntyp_molec(1)))
2609 allocate(rborn_scphoi(ntyp_molec(1)))
2610 allocate(rborn_scphoj(ntyp_molec(1)))
2611 allocate(alphi_scpho(ntyp_molec(1)))
2615 do j=1,ntyp_molec(1) ! without U then we will take T for U
2616 write (*,*) "Im in scpho ", i, " ", j
2617 read(isidep_scpho,*) &
2618 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2619 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2620 write(*,*) "eps",eps_scpho(j)
2621 read(isidep_scpho,*) &
2622 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2623 chis_scpho(j,1),chis_scpho(j,2)
2624 read(isidep_scpho,*) &
2625 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2626 read(isidep_scpho,*) &
2627 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2631 allocate(aa_scpho(ntyp_molec(1)))
2632 allocate(bb_scpho(ntyp_molec(1)))
2633 do j=1,ntyp_molec(1)
2640 sigeps=dsign(1.0D0,epsij)
2642 aa_scpho(j)=epsij*rrij*rrij
2643 bb_scpho(j)=-sigeps*epsij*rrij
2647 read(isidep_peppho,*) &
2648 eps_peppho,sigma_peppho
2649 read(isidep_peppho,*) &
2650 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2651 read(isidep_peppho,*) &
2652 (wqdip_peppho(k),k=1,2)
2660 sigeps=dsign(1.0D0,epsij)
2662 aa_peppho=epsij*rrij*rrij
2663 bb_peppho=-sigeps*epsij*rrij
2666 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2674 ! Define the SC-p interaction constants
2683 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2684 read (itorp_nucl,*) ntortyp_nucl
2685 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2686 !el from energy module---------
2687 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2688 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2690 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2691 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2692 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2693 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2695 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2696 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2697 !el---------------------------
2700 !el--------------------
2701 read (itorp_nucl,*) &
2702 (itortyp_nucl(i),i=1,ntyp_molec(2))
2703 ! print *,itortyp_nucl(:)
2704 !c write (iout,*) 'ntortyp',ntortyp
2707 read (itorp_nucl,*) nterm_nucl(i,j),nlor_nucl(i,j)
2708 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2711 do k=1,nterm_nucl(i,j)
2712 read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2713 v0ij=v0ij+si*v1_nucl(k,i,j)
2716 do k=1,nlor_nucl(i,j)
2717 read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),&
2718 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2719 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2726 !elwrite(iout,*) "parmread kontrol before oldscp"
2728 ! Define the SC-p interaction constants
2732 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2734 ! aad(i,1)=0.3D0*4.0D0**12
2735 ! Following line for constants currently implemented
2736 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2737 aad(i,1)=1.5D0*4.0D0**12
2738 ! aad(i,1)=0.17D0*5.6D0**12
2740 ! "Soft" SC-p repulsion
2742 ! Following line for constants currently implemented
2743 ! aad(i,1)=0.3D0*4.0D0**6
2744 ! "Hard" SC-p repulsion
2745 bad(i,1)=3.0D0*4.0D0**6
2746 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2755 ! 8/9/01 Read the SC-p interaction constants from file
2758 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
2761 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2762 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2763 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2764 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2768 write (iout,*) "Parameters of SC-p interactions:"
2770 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2771 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2775 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2777 do i=1,ntyp_molec(2)
2778 read (iscpp_nucl,*) eps_scp_nucl(i),rscp_nucl(i)
2780 do i=1,ntyp_molec(2)
2781 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2782 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2784 r0pp=1.12246204830937298142*5.16158
2790 ! Define the constants of the disulfide bridge
2794 ! Old arbitrary potential - commented out.
2799 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2800 ! energy surface of diethyl disulfide.
2801 ! A. Liwo and U. Kozlowska, 11/24/03
2812 write (iout,'(/a)') "Disulfide bridge parameters:"
2813 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2814 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2815 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2816 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2820 end subroutine parmread
2822 !-----------------------------------------------------------------------------
2824 !-----------------------------------------------------------------------------
2825 subroutine mygetenv(string,var)
2829 ! This subroutine passes the environmental variables to FORTRAN program.
2830 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
2831 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
2832 ! reads the environmental variables from $HOME/.env
2834 ! Usage: As for the standard FORTRAN GETENV subroutine.
2836 ! Purpose: some versions/installations of MPI do not transfer the environmental
2837 ! variables to slave processors, if these variables are set in the shell script
2838 ! from which mpirun is called.
2847 character*(*) :: string,var
2848 #if defined(MYGETENV) && defined(MPI)
2849 ! include "DIMENSIONS.ZSCOPT"
2851 ! include "COMMON.MPI"
2852 !el character*360 ucase
2854 character(len=360) :: string1(360),karta
2855 character(len=240) :: home
2858 call getenv("HOME",home)
2859 open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
2861 read (99,end=111,err=111,'(a)') karta
2865 call split_string(karta,string1,80,n)
2866 if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
2867 string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
2869 print *,"Processor",me,": ",var(:ilen(var)),&
2870 " assigned to ",string(:ilen(string))
2875 111 print *,"Environment variable ",string(:ilen(string))," not set."
2878 112 print *,"Error opening environment file!"
2880 call getenv(string,var)
2883 end subroutine mygetenv
2884 !-----------------------------------------------------------------------------
2886 !-----------------------------------------------------------------------------
2887 subroutine read_general_data(*)
2889 use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
2890 scelemode,TUBEmode,tor_mode
2892 use energy_data, only:distchainmax,tubeR0,tubecenter
2893 use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
2894 bordtubebot,tubebufthick,buftubebot,buftubetop
2896 ! include "DIMENSIONS"
2897 ! include "DIMENSIONS.ZSCOPT"
2898 ! include "DIMENSIONS.FREE"
2899 ! include "COMMON.TORSION"
2900 ! include "COMMON.INTERACT"
2901 ! include "COMMON.IOUNITS"
2902 ! include "COMMON.TIME1"
2903 ! include "COMMON.PROT"
2904 ! include "COMMON.PROTFILES"
2905 ! include "COMMON.CHAIN"
2906 ! include "COMMON.NAMES"
2907 ! include "COMMON.FFIELD"
2908 ! include "COMMON.ENEPS"
2909 ! include "COMMON.WEIGHTS"
2910 ! include "COMMON.FREE"
2911 ! include "COMMON.CONTROL"
2912 ! include "COMMON.ENERGIES"
2913 character(len=800) :: controlcard
2914 integer :: i,j,k,ii,n_ene_found
2915 integer :: ind,itype1,itype2,itypf,itypsc,itypp
2918 !el character*16 ucase
2919 character(len=16) :: key
2921 call card_concat(controlcard,.true.)
2922 call readi(controlcard,"N_ENE",n_eneW,max_eneW)
2923 if (n_eneW.gt.max_eneW) then
2924 write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
2928 call readi(controlcard,"NPARMSET",nparmset,1)
2929 !elwrite(iout,*)"in read_gen data"
2930 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
2931 call readi(controlcard,"IPARMPRINT",iparmprint,1)
2932 write (iout,*) "PARMPRINT",iparmprint
2933 if (nparmset.gt.max_parm) then
2934 write (iout,*) "Error: parameter out of range: NPARMSET",&
2938 !elwrite(iout,*)"in read_gen data"
2939 call readi(controlcard,"MAXIT",maxit,5000)
2940 call reada(controlcard,"FIMIN",fimin,1.0d-3)
2941 call readi(controlcard,"ENSEMBLES",ensembles,0)
2942 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
2943 write (iout,*) "Number of energy parameter sets",nparmset
2944 allocate(isampl(nparmset))
2945 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
2946 write (iout,*) "MaxSlice",MaxSlice
2947 call readi(controlcard,"NSLICE",nslice,1)
2948 !elwrite(iout,*)"in read_gen data"
2950 if (nslice.gt.MaxSlice) then
2951 write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
2955 write (iout,*) "Frequency of storing conformations",&
2956 (isampl(i),i=1,nparmset)
2957 write (iout,*) "Maxit",maxit," Fimin",fimin
2958 call readi(controlcard,"NQ",nQ,1)
2959 if (nQ.gt.MaxQ) then
2960 write (iout,*) "Error: parameter out of range: NQ",nq,&
2965 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
2966 call reada(controlcard,"DELTA",delta,1.0d-2)
2967 call readi(controlcard,"EINICHECK",einicheck,2)
2968 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
2969 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
2970 call readi(controlcard,"RESCALE",rescale_modeW,1)
2971 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2972 call reada(controlcard,'BOXY',boxysize,100.0d0)
2973 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2974 call readi(controlcard,"SCELEMODE",scelemode,0)
2975 print *,"SCELE",scelemode
2976 call readi(controlcard,'TORMODE',tor_mode,0)
2977 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2978 write(iout,*) "torsional and valence angle mode",tor_mode
2980 call readi(controlcard,'TUBEMOD',tubemode,0)
2982 if (TUBEmode.gt.0) then
2983 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
2984 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
2985 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
2986 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
2987 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
2988 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
2989 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
2990 buftubebot=bordtubebot+tubebufthick
2991 buftubetop=bordtubetop-tubebufthick
2993 ions=index(controlcard,"IONS").gt.0
2994 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2995 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2996 write(iout,*) "R_CUT_ELE=",r_cut_ele
2997 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
2998 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
2999 call readi(controlcard,'SYM',symetr,1)
3000 write (iout,*) "DISTCHAINMAX",distchainmax
3001 write (iout,*) "delta",delta
3002 write (iout,*) "einicheck",einicheck
3003 write (iout,*) "rescale_mode",rescale_modeW
3005 bxfile=index(controlcard,"BXFILE").gt.0
3006 cxfile=index(controlcard,"CXFILE").gt.0
3007 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
3009 histfile=index(controlcard,"HISTFILE").gt.0
3010 histout=index(controlcard,"HISTOUT").gt.0
3011 entfile=index(controlcard,"ENTFILE").gt.0
3012 zscfile=index(controlcard,"ZSCFILE").gt.0
3013 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
3014 write (iout,*) "with_dihed_constr ",with_dihed_constr
3015 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3017 end subroutine read_general_data
3018 !------------------------------------------------------------------------------
3019 subroutine read_efree(*)
3021 ! Read molecular data
3024 ! include 'DIMENSIONS'
3025 ! include 'DIMENSIONS.ZSCOPT'
3026 ! include 'DIMENSIONS.COMPAR'
3027 ! include 'DIMENSIONS.FREE'
3028 ! include 'COMMON.IOUNITS'
3029 ! include 'COMMON.TIME1'
3030 ! include 'COMMON.SBRIDGE'
3031 ! include 'COMMON.CONTROL'
3032 ! include 'COMMON.CHAIN'
3033 ! include 'COMMON.HEADER'
3034 ! include 'COMMON.GEO'
3035 ! include 'COMMON.FREE'
3036 character(len=320) :: controlcard !,ucase
3037 integer :: iparm,ib,i,j,npars
3047 ! call alloc_wham_arrays
3048 ! allocate(nT_h(nParmSet))
3049 ! allocate(replica(nParmSet))
3050 ! allocate(umbrella(nParmSet))
3051 ! allocate(read_iset(nParmSet))
3052 ! allocate(nT_h(nParmSet))
3056 call card_concat(controlcard,.true.)
3057 call readi(controlcard,'NT',nT_h(iparm),1)
3058 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
3060 if (nT_h(iparm).gt.MaxT_h) then
3061 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
3065 replica(iparm)=index(controlcard,"REPLICA").gt.0
3066 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
3067 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
3068 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
3069 replica(iparm)," umbrella ",umbrella(iparm),&
3070 " read_iset",read_iset(iparm)
3073 call card_concat(controlcard,.true.)
3074 call readi(controlcard,'NR',nR(ib,iparm),1)
3075 if (umbrella(iparm)) then
3078 nRR(ib,iparm)=nR(ib,iparm)
3080 if (nR(ib,iparm).gt.MaxR) then
3081 write (iout,*) "Error: parameter out of range: NR",&
3085 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
3086 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
3087 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
3090 call card_concat(controlcard,.true.)
3091 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
3093 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
3098 write (iout,*) "ib",ib," beta_h",&
3099 1.0d0/(0.001987*beta_h(ib,iparm))
3100 write (iout,*) "nR",nR(ib,iparm)
3101 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
3103 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
3104 "q0",(q0(j,i,ib,iparm),j=1,nQ)
3116 nR(ib,iparm)=nR(ib,1)
3117 if (umbrella(iparm)) then
3120 nRR(ib,iparm)=nR(ib,1)
3122 beta_h(ib,iparm)=beta_h(ib,1)
3124 f(i,ib,iparm)=f(i,ib,1)
3126 KH(j,i,ib,iparm)=KH(j,i,ib,1)
3127 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
3130 replica(iparm)=replica(1)
3131 umbrella(iparm)=umbrella(1)
3132 read_iset(iparm)=read_iset(1)
3139 end subroutine read_efree
3140 !-----------------------------------------------------------------------------
3141 subroutine read_protein_data(*)
3143 ! include "DIMENSIONS"
3144 ! include "DIMENSIONS.ZSCOPT"
3145 ! include "DIMENSIONS.FREE"
3149 integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
3150 ! include "COMMON.MPI"
3152 ! include "COMMON.CHAIN"
3153 ! include "COMMON.IOUNITS"
3154 ! include "COMMON.PROT"
3155 ! include "COMMON.PROTFILES"
3156 ! include "COMMON.NAMES"
3157 ! include "COMMON.FREE"
3158 ! include "COMMON.OBCINKA"
3159 character(len=64) :: nazwa
3160 character(len=16000) :: controlcard
3161 integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
3162 !el external ilen,iroof
3171 ! Read names of files with conformation data.
3172 if (replica(iparm)) then
3178 do ii=1,nRR(ib,iparm)
3179 write (iout,*) "Parameter set",iparm," temperature",ib,&
3182 call card_concat(controlcard,.true.)
3183 write (iout,*) controlcard(:ilen(controlcard))
3184 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
3185 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
3186 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
3187 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
3188 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
3189 maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
3190 call reada(controlcard,"TIME_START",&
3191 time_start_collect(ii,ib,iparm),0.0d0)
3192 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
3194 write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
3195 " rec_end",rec_end(ii,ib,iparm)
3196 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
3197 " time_end",time_end_collect(ii,ib,iparm)
3199 if (replica(iparm)) then
3200 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
3201 write (iout,*) "Number of trajectories",totraj(ii,iparm)
3204 if (nfile_bin(ii,ib,iparm).lt.2 &
3205 .and. nfile_asc(ii,ib,iparm).eq.0 &
3206 .and. nfile_cx(ii,ib,iparm).eq.0) then
3207 write (iout,*) "Error - no action specified!"
3210 if (nfile_bin(ii,ib,iparm).gt.0) then
3211 call card_concat(controlcard,.false.)
3212 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
3213 maxfile_prot,nfile_bin(ii,ib,iparm))
3215 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
3216 write(iout,*) (protfiles(i,1,ii,ib,iparm),&
3217 i=1,nfile_bin(ii,ib,iparm))
3220 if (nfile_asc(ii,ib,iparm).gt.0) then
3221 call card_concat(controlcard,.false.)
3222 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3223 maxfile_prot,nfile_asc(ii,ib,iparm))
3225 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
3226 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3227 i=1,nfile_asc(ii,ib,iparm))
3229 else if (nfile_cx(ii,ib,iparm).gt.0) then
3230 call card_concat(controlcard,.false.)
3231 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3232 maxfile_prot,nfile_cx(ii,ib,iparm))
3234 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
3235 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3236 i=1,nfile_cx(ii,ib,iparm))
3245 end subroutine read_protein_data
3246 !-------------------------------------------------------------------------------
3247 subroutine readsss(rekord,lancuch,wartosc,default)
3249 character*(*) :: rekord,lancuch,wartosc,default
3250 character(len=80) :: aux
3251 integer :: lenlan,lenrec,iread,ireade
3255 lenlan=ilen(lancuch)
3257 iread=index(rekord,lancuch(:lenlan)//"=")
3258 ! print *,"rekord",rekord," lancuch",lancuch
3259 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3260 if (iread.eq.0) then
3264 iread=iread+lenlan+1
3265 ! print *,"iread",iread
3266 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3267 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3269 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3271 ! print *,"iread",iread
3272 if (iread.gt.lenrec) then
3277 ! print *,"ireade",ireade
3278 do while (ireade.lt.lenrec .and. &
3279 .not.iblnk(rekord(ireade:ireade)))
3282 wartosc=rekord(iread:ireade)
3284 end subroutine readsss
3285 !----------------------------------------------------------------------------
3286 subroutine multreads(rekord,lancuch,tablica,dim,default)
3289 character*(*) rekord,lancuch,tablica(dim),default
3290 character(len=80) :: aux
3291 integer :: lenlan,lenrec,iread,ireade
3298 lenlan=ilen(lancuch)
3300 iread=index(rekord,lancuch(:lenlan)//"=")
3301 ! print *,"rekord",rekord," lancuch",lancuch
3302 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3303 if (iread.eq.0) return
3304 iread=iread+lenlan+1
3306 ! print *,"iread",iread
3307 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3308 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3310 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3312 ! print *,"iread",iread
3313 if (iread.gt.lenrec) return
3315 ! print *,"ireade",ireade
3316 do while (ireade.lt.lenrec .and. &
3317 .not.iblnk(rekord(ireade:ireade)))
3320 tablica(i)=rekord(iread:ireade)
3323 end subroutine multreads
3324 !----------------------------------------------------------------------------
3325 subroutine split_string(rekord,tablica,dim,nsub)
3327 integer :: dim,nsub,i,ii,ll,kk
3328 character*(*) tablica(dim)
3329 character*(*) rekord
3339 ! Find the start of term name
3341 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
3344 ! Parse the name into TABLICA(i) until blank found
3345 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
3347 tablica(i)(kk:kk)=rekord(ii:ii)
3350 if (kk.gt.0) nsub=nsub+1
3351 if (ii.gt.ll) return
3354 end subroutine split_string
3355 !--------------------------------------------------------------------------------
3357 !--------------------------------------------------------------------------------
3358 subroutine read_compar
3360 ! Read molecular data
3362 use conform_compar, only:alloc_compar_arrays
3363 use control_data, only:pdbref
3364 use geometry_data, only:deg2rad,rad2deg
3366 ! include 'DIMENSIONS'
3367 ! include 'DIMENSIONS.ZSCOPT'
3368 ! include 'DIMENSIONS.COMPAR'
3369 ! include 'DIMENSIONS.FREE'
3370 ! include 'COMMON.IOUNITS'
3371 ! include 'COMMON.TIME1'
3372 ! include 'COMMON.SBRIDGE'
3373 ! include 'COMMON.CONTROL'
3374 ! include 'COMMON.COMPAR'
3375 ! include 'COMMON.CHAIN'
3376 ! include 'COMMON.HEADER'
3377 ! include 'COMMON.GEO'
3378 ! include 'COMMON.FREE'
3379 character(len=320) :: controlcard !,ucase
3380 character(len=64) :: wfile
3384 !elwrite(iout,*)"jestesmy w read_compar"
3385 call card_concat(controlcard,.true.)
3386 pdbref=(index(controlcard,'PDBREF').gt.0)
3387 call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
3388 call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
3389 call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
3390 call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
3391 verbose = index(controlcard,"VERBOSE").gt.0
3392 lgrp=index(controlcard,"STATIN").gt.0
3393 lgrp_out=index(controlcard,"STATOUT").gt.0
3394 merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
3395 binary = index(controlcard,"BINARY").gt.0
3396 rmscut_base_up=rmscut_base_up/50
3397 rmscut_base_low=rmscut_base_low/50
3398 call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
3399 call readi(controlcard,'NLEVEL',nlevel,1)
3400 if (nlevel.lt.0) then
3402 call alloc_compar_arrays(maxfrag,1)
3405 allocate(nfrag(nlevel))
3407 ! Read the data pertaining to elementary fragments (level 1)
3408 call readi(controlcard,'NFRAG',nfrag(1),0)
3409 write(iout,*)"nfrag(1)",nfrag(1)
3410 call alloc_compar_arrays(nfrag(1),nlevel)
3412 call card_concat(controlcard,.true.)
3413 write (iout,*) controlcard(:ilen(controlcard))
3414 call readi(controlcard,'NPIECE',npiece(j,1),0)
3415 call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
3416 call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
3417 call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
3418 call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
3419 call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
3420 call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
3421 call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
3422 call readi(controlcard,'RMS',irms(j,1),0)
3423 call readi(controlcard,'LOCAL',iloc(j),1)
3424 call readi(controlcard,'ELCONT',ielecont(j,1),1)
3425 if (ielecont(j,1).eq.0) then
3426 call readi(controlcard,'SCCONT',isccont(j,1),1)
3428 ang_cut(j)=ang_cut(j)*deg2rad
3429 ang_cut1(j)=ang_cut1(j)*deg2rad
3431 call card_concat(controlcard,.true.)
3432 call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
3433 call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
3435 write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
3436 (ifrag(1,k,j),ifrag(2,k,j),&
3437 k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
3438 " ang_cut1",ang_cut1(j)*rad2deg
3439 write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
3440 write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
3441 write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
3442 " ilocal",iloc(j)," isccont",isccont(j,1)
3444 ! Read data pertaning to higher levels
3446 call card_concat(controlcard,.true.)
3447 call readi(controlcard,'NFRAG',NFRAG(i),0)
3448 write (iout,*) "i",i," nfrag",nfrag(i)
3450 call card_concat(controlcard,.true.)
3452 call readi(controlcard,'ELCONT',ielecont(j,i),0)
3453 if (ielecont(j,i).eq.0) then
3454 call readi(controlcard,'SCCONT',isccont(j,i),1)
3456 call readi(controlcard,'RMS',irms(j,i),0)
3462 call readi(controlcard,'NPIECE',npiece(j,i),0)
3463 call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
3464 call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
3465 call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
3467 call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
3468 call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
3469 write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
3470 n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
3471 " isccont",isccont(j,i)," irms",irms(j,i)
3472 write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
3473 write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
3474 write(iout,*)"nc_frac",nc_fragm(j,i),&
3475 " nc_req",nc_req_setf(j,i)
3478 if (binary) write (iout,*) "Classes written in binary format."
3481 call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
3482 call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
3483 call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
3484 call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
3485 call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
3486 call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
3487 call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
3488 call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
3489 call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
3490 call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
3491 call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
3492 call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
3493 call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
3494 call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
3495 call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
3496 call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
3497 call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
3498 call readi(controlcard,'RMS_SINGLE',irms_single,0)
3499 call readi(controlcard,'CONT_SINGLE',icont_single,1)
3500 call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
3501 call readi(controlcard,'RMS_PAIR',irms_pair,0)
3502 call readi(controlcard,'CONT_PAIR',icont_pair,1)
3503 call readi(controlcard,'SPLIT_BET',isplit_bet,0)
3504 angcut_hel=angcut_hel*deg2rad
3505 angcut1_hel=angcut1_hel*deg2rad
3506 angcut_bet=angcut_bet*deg2rad
3507 angcut1_bet=angcut1_bet*deg2rad
3508 angcut_strand=angcut_strand*deg2rad
3509 angcut1_strand=angcut1_strand*deg2rad
3510 write (iout,*) "Automatic detection of structural elements"
3511 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
3512 ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
3513 ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
3514 ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
3515 ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
3516 ' SPLIT_BET',isplit_bet
3517 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
3518 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
3519 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
3520 ' MAXANG_HEL',angcut1_hel*rad2deg
3521 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
3522 ' MAXANG_BET',angcut1_bet*rad2deg
3523 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
3524 ' MAXANG_STRAND',angcut1_strand*rad2deg
3525 write (iout,*) 'FRAC_MIN',frac_min_set
3527 end subroutine read_compar
3528 !--------------------------------------------------------------------------------
3530 !--------------------------------------------------------------------------------
3531 subroutine read_ref_structure(*)
3533 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
3536 use control_data, only:pdbref
3537 use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
3538 nstart_sup,nstart_seq,nperm,nres0
3539 use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
3540 use compare, only:seq_comp !,contact,elecont
3541 use geometry, only:chainbuild,dist
3542 use io_config, only:readpdb
3544 use conform_compar, only:contact,elecont
3546 ! include 'DIMENSIONS'
3547 ! include 'DIMENSIONS.ZSCOPT'
3548 ! include 'DIMENSIONS.COMPAR'
3549 ! include 'COMMON.IOUNITS'
3550 ! include 'COMMON.GEO'
3551 ! include 'COMMON.VAR'
3552 ! include 'COMMON.INTERACT'
3553 ! include 'COMMON.LOCAL'
3554 ! include 'COMMON.NAMES'
3555 ! include 'COMMON.CHAIN'
3556 ! include 'COMMON.FFIELD'
3557 ! include 'COMMON.SBRIDGE'
3558 ! include 'COMMON.HEADER'
3559 ! include 'COMMON.CONTROL'
3560 ! include 'COMMON.CONTACTS1'
3561 ! include 'COMMON.PEPTCONT'
3562 ! include 'COMMON.TIME1'
3563 ! include 'COMMON.COMPAR'
3564 character(len=4) :: sequence(nres)
3566 !el real(kind=8) :: x(maxvar)
3567 integer :: itype_pdb(nres,5)
3568 !el logical seq_comp
3569 integer :: i,j,k,nres_pdb,iaux,mnum
3570 real(kind=8) :: ddsc !el,dist
3571 integer :: kkk !,ilen
3575 write (iout,*) "pdbref",pdbref
3577 read(inp,'(a)') pdbfile
3578 write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
3579 pdbfile(:ilen(pdbfile))
3580 open(ipdbin,file=pdbfile,status='old',err=33)
3582 33 write (iout,'(a)') 'Error opening PDB file.'
3587 itype_pdb(i,mnum)=itype(i,mnum)
3593 iaux=itype_pdb(i,mnum)
3594 itype_pdb(i,mnum)=itype(i,mnum)
3602 if (nsup.le.(nct-nnt+1)) then
3603 do i=0,nct-nnt+1-nsup
3604 if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
3606 do j=nnt+nsup-1,nnt,-1
3608 cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
3611 do j=nnt+nsup-1,nnt,-1
3613 cref(k,j+i,kkk)=cref(k,j,kkk)
3615 write(iout,*) "J",j,"J+I",j+i
3616 phi_ref(j+i)=phi_ref(j)
3617 theta_ref(j+i)=theta_ref(j)
3618 alph_ref(j+i)=alph_ref(j)
3619 omeg_ref(j+i)=omeg_ref(j)
3623 write (iout,'(i5,3f10.5,5x,3f10.5)') &
3624 j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
3632 write (iout,'(a)') &
3633 'Error - sequences to be superposed do not match.'
3636 do i=0,nsup-(nct-nnt+1)
3637 if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
3640 nstart_sup=nstart_sup+i
3645 write (iout,'(a)') &
3646 'Error - sequences to be superposed do not match.'
3650 write (iout,'(a,i5)') &
3651 'Experimental structure begins at residue',nstart_seq
3653 call read_angles(inp,*38)
3655 38 write (iout,'(a)') 'Error reading reference structure.'
3664 cref(j,i,kkk)=c(j,i)
3668 nend_sup=nstart_sup+nsup-1
3671 c(j,i)=cref(j,i,kkk)
3677 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
3679 if (itype(i,mnum).ne.10) then
3680 ddsc = dist(i,nres+i)
3682 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
3686 dc_norm(j,nres+i)=0.0d0
3689 ! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
3690 ! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
3691 ! dc_norm(3,nres+i)**2
3693 dc(j,i)=c(j,i+1)-c(j,i)
3697 dc_norm(j,i)=dc(j,i)/ddsc
3700 ! print *,"Calling contact"
3701 call contact(.true.,ncont_ref,icont_ref(1,1),&
3702 nstart_sup,nend_sup)
3703 ! print *,"Calling elecont"
3704 call elecont(.true.,ncont_pept_ref,&
3705 icont_pept_ref(1,1),&
3706 nstart_sup,nend_sup)
3707 write (iout,'(a,i3,a,i3,a,i3,a)') &
3708 'Number of residues to be superposed:',nsup,&
3709 ' (from residue',nstart_sup,' to residue',&
3712 end subroutine read_ref_structure
3713 !--------------------------------------------------------------------------------
3715 !--------------------------------------------------------------------------------
3716 subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
3718 use geometry_data, only:nres,c
3719 use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
3720 ! implicit real*8 (a-h,o-z)
3721 ! include 'DIMENSIONS'
3722 ! include 'DIMENSIONS.ZSCOPT'
3723 ! include 'COMMON.CHAIN'
3724 ! include 'COMMON.INTERACT'
3725 ! include 'COMMON.NAMES'
3726 ! include 'COMMON.IOUNITS'
3727 ! include 'COMMON.HEADER'
3728 ! include 'COMMON.SBRIDGE'
3729 character(len=50) :: tytul
3730 character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',&
3731 'D','E','F','G','H','I','J','K','L','M','N','O',&
3732 'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
3733 integer,dimension(nres) :: ica !(maxres)
3734 real(kind=8) :: temp,efree,etot,entropy,rmsdev
3735 integer :: ii,i,j,iti,ires,iatom,ichain,mnum
3736 write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
3738 write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
3740 write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
3748 if (iti.eq.ntyp1) then
3749 if (itype(i-1,molnum(i-1)).eq.ntyp1) then
3752 write (ipdb,'(a)') 'TER'
3758 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
3762 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
3763 ires,(c(j,nres+i),j=1,3)
3767 write (ipdb,'(a)') 'TER'
3770 if (itype(i,mnum).eq.ntyp1) cycle
3771 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3772 write (ipdb,30) ica(i),ica(i+1)
3773 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3774 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
3775 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
3776 write (ipdb,30) ica(i),ica(i)+1
3779 if (itype(nct,molnum(nct)).ne.10) then
3780 write (ipdb,30) ica(nct),ica(nct)+1
3783 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
3785 write (ipdb,'(a6)') 'ENDMDL'
3786 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3787 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3788 30 FORMAT ('CONECT',8I5)
3790 end subroutine pdboutW
3792 !------------------------------------------------------------------------------
3794 !-----------------------------------------------------------------------------
3795 !-----------------------------------------------------------------------------