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')
101 call mygetenv('IONPAR_NUCL',ionnuclname)
102 open (iionnucl,file=ionnuclname,status='old')
107 ! 8/9/01 In the newest version SCp interaction constants are read from a file
108 ! Use -DOLDSCP to use hard-coded constants instead.
110 call mygetenv('SCPPAR',scpname)
111 open (iscpp,file=scpname,status='old')
114 if (MyID.eq.BossID) then
115 MyRank = MyID/fgProcs
118 print *,'OpenUnits: processor',MyRank
119 call numstr(MyRank,liczba)
120 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//liczba
122 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
124 open(iout,file=outname,status='unknown')
125 write (iout,'(80(1h-))')
126 write (iout,'(30x,a)') "FILE ASSIGNMENT"
127 write (iout,'(80(1h-))')
128 write (iout,*) "Input file : ",&
129 prefix(:ilen(prefix))//'.inp'
130 write (iout,*) "Output file : ",&
131 outname(:ilen(outname))
133 write (iout,*) "Sidechain potential file : ",&
134 sidename(:ilen(sidename))
136 write (iout,*) "SCp potential file : ",&
137 scpname(:ilen(scpname))
139 write (iout,*) "Electrostatic potential file : ",&
140 elename(:ilen(elename))
141 write (iout,*) "Cumulant coefficient file : ",&
142 fouriername(:ilen(fouriername))
143 write (iout,*) "Torsional parameter file : ",&
144 torname(:ilen(torname))
145 write (iout,*) "Double torsional parameter file : ",&
146 tordname(:ilen(tordname))
147 write (iout,*) "Backbone-rotamer parameter file : ",&
148 sccorname(:ilen(sccorname))
149 write (iout,*) "Bond & inertia constant file : ",&
150 bondname(:ilen(bondname))
151 write (iout,*) "Bending parameter file : ",&
152 thetname(:ilen(thetname))
153 write (iout,*) "Rotamer parameter file : ",&
154 rotname(:ilen(rotname))
155 write (iout,'(80(1h-))')
158 end subroutine openunits
159 !-----------------------------------------------------------------------------
161 !-----------------------------------------------------------------------------
162 subroutine molread(*)
164 ! Read molecular data.
167 use geometry_data, only:nres,deg2rad,c,dc,nres_molec
168 use control_data, only:iscode,dyn_ss
169 use io_base, only:rescode
170 use control, only:setup_var,init_int_table,hpb_partition
171 use geometry, only:alloc_geo_arrays
172 use energy, only:alloc_ener_arrays
173 ! implicit real*8 (a-h,o-z)
174 ! include 'DIMENSIONS'
175 ! include 'DIMENSIONS.ZSCOPT'
176 ! include 'COMMON.IOUNITS'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.VAR'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.LOCAL'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.CHAIN'
183 ! include 'COMMON.FFIELD'
184 ! include 'COMMON.SBRIDGE'
185 ! include 'COMMON.TORCNSTR'
186 ! include 'COMMON.CONTROL'
187 character(len=4),dimension(:),allocatable :: sequence !(nres)
188 !el integer :: rescode
189 !el real(kind=8) :: x(maxvar)
190 character(len=320) :: controlcard !,ucase
191 integer,dimension(nres,5) :: itype_pdb !(maxres)
192 integer :: i,j,i1,i2,it1,it2,mnum
193 real(kind=8) :: scalscp
194 !el logical :: seq_comp
195 write(iout,*) "START MOLREAD"
199 print *,"nres_molec, initial",nres_molec(i)
202 call card_concat(controlcard,.true.)
203 call reada(controlcard,'SCAL14',scal14,0.4d0)
204 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
205 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
206 call reada(controlcard,'TEMP0',temp0,300.0d0) !el
207 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
208 r0_corr=cutoff_corr-delt_corr
209 call readi(controlcard,"NRES",nres_molec(1),0)
210 call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
211 call readi(controlcard,"NRES_CAT",nres_molec(5),0)
214 nres=nres_molec(i)+nres
216 ! allocate(sequence(nres+1))
217 !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
218 call alloc_geo_arrays
219 call alloc_ener_arrays
220 ! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
221 !----------------------------
222 allocate(c(3,2*nres+2))
223 allocate(dc(3,0:2*nres+2))
224 allocate(itype(nres+2,5))
225 allocate(itel(nres+2))
226 if (.not. allocated(molnum)) allocate(molnum(nres+2))
242 !--------------------------
244 allocate(sequence(nres+1))
245 iscode=index(controlcard,"ONE_LETTER")
247 write (iout,*) "Error: no residues in molecule"
250 if (nres.gt.maxres) then
251 write (iout,*) "Error: too many residues",nres,maxres
253 write(iout,*) 'nres=',nres
255 ! Read sequence of the protein
256 if (iscode.gt.0) then
257 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
259 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
261 ! Convert sequence to numeric code
265 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
267 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
270 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
272 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
275 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
278 write (iout,*) "Numeric code:"
279 write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres)
283 if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
285 if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
289 else if (iabs(itype(i+1,mnum)).ne.20) then
291 else if (iabs(itype(i,mnum)).ne.20) then
298 write (iout,*) "ITEL"
301 write (iout,*) i,itype(i,mnum),itel(i)
306 if (with_dihed_constr) then
308 read (inp,*) ndih_constr
309 if (ndih_constr.gt.0) then
311 write (iout,*) 'FTORS',ftors
312 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
314 'There are',ndih_constr,' constraints on phi angles.'
316 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
319 phi0(i)=deg2rad*phi0(i)
320 drange(i)=deg2rad*drange(i)
328 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
329 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
330 write(iout,*) 'NNT=',NNT,' NCT=',NCT
334 write (iout,'(/a,i3,a)') 'The chain contains',ns,&
335 ' disulfide-bridging cysteines.'
336 write (iout,'(20i4)') (iss(i),i=1,ns)
337 write (iout,'(/a/)') 'Pre-formed links are:'
344 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
345 restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',&
346 dhpb(i),ebr,forcon(i)
349 if (ns.gt.0.and.dyn_ss) then
353 forcon(i-nss)=forcon(i)
360 dyn_ss_mask(iss(i))=.true.
365 end subroutine molread
366 !-----------------------------------------------------------------------------
368 !-----------------------------------------------------------------------------
369 subroutine parmread(iparm,*)
374 ! Read the parameters of the probability distributions of the virtual-bond
375 ! valence angles and the side chains and energy parameters.
381 use control_data, only: maxterm,maxlor,maxterm_sccor,& !maxtor
382 maxtermd_1,maxtermd_2,tor_mode,scelemode !,maxthetyp,maxthetyp1
386 use io_config, only: printmat
387 use control, only: getenv_loc
394 ! implicit real*8 (a-h,o-z)
395 ! include 'DIMENSIONS'
396 ! include 'DIMENSIONS.ZSCOPT'
397 ! include 'DIMENSIONS.FREE'
398 ! include 'COMMON.IOUNITS'
399 ! include 'COMMON.CHAIN'
400 ! include 'COMMON.INTERACT'
401 ! include 'COMMON.GEO'
402 ! include 'COMMON.LOCAL'
403 ! include 'COMMON.TORSION'
404 ! include 'COMMON.FFIELD'
405 ! include 'COMMON.NAMES'
406 ! include 'COMMON.SBRIDGE'
407 ! include 'COMMON.WEIGHTS'
408 ! include 'COMMON.ENEPS'
409 ! include 'COMMON.SCCOR'
410 ! include 'COMMON.SCROT'
411 ! include 'COMMON.FREE'
412 character(len=1) :: t1,t2,t3
413 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
414 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
415 logical :: lprint,SPLIT_FOURIERTOR
416 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
417 character(len=800) :: controlcard
418 character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
419 tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
423 character(len=16) :: key
424 integer :: iparm,nkcctyp
425 !el real(kind=8) :: ip,mp
426 real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
427 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm, &
428 epspeptube,epssctube,sigmapeptube, &
431 real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
433 integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj
434 integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
440 ! Set LPRINT=.TRUE. for debugging
441 dwa16=2.0d0**(1.0d0/6.0d0)
444 ! Assign virtual-bond length
447 vblinv2=vblinv*vblinv
450 call card_concat(controlcard,.true.)
453 call reada(controlcard,"D0CM",d0cm,3.78d0)
454 call reada(controlcard,"AKCM",akcm,15.1d0)
455 call reada(controlcard,"AKTH",akth,11.0d0)
456 call reada(controlcard,"AKCT",akct,12.0d0)
457 call reada(controlcard,"V1SS",v1ss,-1.08d0)
458 call reada(controlcard,"V2SS",v2ss,7.61d0)
459 call reada(controlcard,"V3SS",v3ss,13.7d0)
460 call reada(controlcard,"EBR",ebr,-5.50D0)
461 call reada(controlcard,"ATRISS",atriss,0.301D0)
462 call reada(controlcard,"BTRISS",btriss,0.021D0)
463 call reada(controlcard,"CTRISS",ctriss,1.001D0)
464 call reada(controlcard,"DTRISS",dtriss,1.001D0)
465 call reada(controlcard,"SSSCALE",ssscale,1.0D0)
466 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
468 allocate(ww(max_eneW))
470 key = wname(i)(:ilen(wname(i)))
471 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
474 write (iout,*) "iparm",iparm," myparm",myparm
475 ! If reading not own parameters, skip assignment
477 if (iparm.eq.myparm .or. .not.separate_parset) then
480 ! Setup weights for UNRES
520 ! print *,"KURWA",ww(48)
521 ! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO "
522 ! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",&
523 ! "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",&
524 ! "WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",&
525 ! "WCATPROT ","WCATCAT
526 ! +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
527 ! +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
528 ! +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
529 ! +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
534 allocate(weights(n_ene))
549 weights(15)=wstrain !0
552 weights(18)=0 !scal14 !
556 weights(26)= wvdwpp_nucl
563 weights(32) =wbond_nucl
564 weights(33) =wang_nucl
566 weights(35) =wtor_nucl
567 weights(36) =wtor_d_nucl
568 weights(37) =wcorr_nucl
569 weights(38) =wcorr3_nucl
571 weights(42) =wcatprot
573 weights(47)= wpepbase
576 weights(50) =wcatnucl
579 call card_concat(controlcard,.false.)
581 ! Return if not own parameters
583 if (iparm.ne.myparm .and. separate_parset) return
585 call reads(controlcard,"BONDPAR",bondname_t,bondname)
586 open (ibond,file=bondname_t,status='old')
588 call reads(controlcard,"THETPAR",thetname_t,thetname)
589 open (ithep,file=thetname_t,status='old')
591 call reads(controlcard,"ROTPAR",rotname_t,rotname)
592 open (irotam,file=rotname_t,status='old')
594 call reads(controlcard,"TORPAR",torname_t,torname)
595 open (itorp,file=torname_t,status='old')
597 call reads(controlcard,"TORDPAR",tordname_t,tordname)
598 open (itordp,file=tordname_t,status='old')
600 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
601 open (isccor,file=sccorname_t,status='old')
603 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
604 open (ifourier,file=fouriername_t,status='old')
606 call reads(controlcard,"ELEPAR",elename_t,elename)
607 open (ielep,file=elename_t,status='old')
609 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
610 open (isidep,file=sidename_t,status='old')
612 call reads(controlcard,"SCPPAR",scpname_t,scpname)
613 open (iscpp,file=scpname_t,status='old')
615 call getenv_loc('IONPAR',ionname)
616 open (iion,file=ionname,status='old')
618 write (iout,*) "Parameter set:",iparm
619 write (iout,*) "Energy-term weights:"
621 write (iout,'(i3,a16,f10.5)') i,wname(i),ww(i)
623 write (iout,*) "Sidechain potential file : ",&
624 sidename_t(:ilen(sidename_t))
626 write (iout,*) "SCp potential file : ",&
627 scpname_t(:ilen(scpname_t))
629 write (iout,*) "Electrostatic potential file : ",&
630 elename_t(:ilen(elename_t))
631 write (iout,*) "Cumulant coefficient file : ",&
632 fouriername_t(:ilen(fouriername_t))
633 write (iout,*) "Torsional parameter file : ",&
634 torname_t(:ilen(torname_t))
635 write (iout,*) "Double torsional parameter file : ",&
636 tordname_t(:ilen(tordname_t))
637 write (iout,*) "Backbone-rotamer parameter file : ",&
638 sccorname(:ilen(sccorname))
639 write (iout,*) "Bond & inertia constant file : ",&
640 bondname_t(:ilen(bondname_t))
641 write (iout,*) "Bending parameter file : ",&
642 thetname_t(:ilen(thetname_t))
643 write (iout,*) "Rotamer parameter file : ",&
644 rotname_t(:ilen(rotname_t))
647 ! Read the virtual-bond parameters, masses, and moments of inertia
648 ! and Stokes' radii of the peptide group and side chains
650 allocate(dsc(ntyp1)) !(ntyp1)
651 allocate(dsc_inv(ntyp1)) !(ntyp1)
652 allocate(nbondterm(ntyp)) !(ntyp)
653 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
654 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
655 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
656 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
657 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
659 !el allocate(msc(ntyp+1)) !(ntyp+1)
660 !el allocate(isc(ntyp+1)) !(ntyp+1)
661 !el allocate(restok(ntyp+1)) !(ntyp+1)
662 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
665 read (ibond,*) vbldp0,akp
668 read (ibond,*) vbldsc0(1,i),aksc(1,i)
669 dsc(i) = vbldsc0(1,i)
673 dsc_inv(i)=1.0D0/dsc(i)
677 read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
679 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
681 dsc(i) = vbldsc0(1,i)
685 dsc_inv(i)=1.0D0/dsc(i)
690 write(iout,'(/a/)')"Force constants virtual bonds:"
691 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
693 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
695 write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
696 vbldsc0(1,i),aksc(1,i),abond0(1,i)
698 write (iout,'(13x,3f10.5)') &
699 vbldsc0(j,i),aksc(j,i),abond0(j,i)
703 if (.not. allocated(msc)) allocate(msc(ntyp1,5))
704 if (.not. allocated(restok)) allocate(restok(ntyp1,5))
705 if (oldion.eq.1) then
708 read(iion,*) msc(i,5),restok(i,5)
709 print *,msc(i,5),restok(i,5)
713 read (iion,*) ncatprotparm
714 allocate(catprm(ncatprotparm,4))
716 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
720 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
724 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
727 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
728 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
732 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
733 restyp(i,5),"-",restyp(j,2)
737 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
740 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2)
741 ! dsc(i) = vbldsc0_nucl(1,i)
745 ! dsc_inv(i)=1.0D0/dsc(i)
750 !----------------------------------------------------
751 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
752 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
753 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
754 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
755 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
756 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
762 athet(j,i,ichir1,ichir2)=0.0D0
763 bthet(j,i,ichir1,ichir2)=0.0D0
777 !elwrite(iout,*) "parmread kontrol"
781 ! Read the parameters of the probability distribution/energy expression
782 ! of the virtual-bond valence angles theta
785 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
786 (bthet(j,i,1,1),j=1,2)
787 read (ithep,*) (polthet(j,i),j=0,3)
788 !elwrite(iout,*) "parmread kontrol in cryst_theta"
789 read (ithep,*) (gthet(j,i),j=1,3)
790 !elwrite(iout,*) "parmread kontrol in cryst_theta"
791 read (ithep,*) theta0(i),sig0(i),sigc0(i)
793 !elwrite(iout,*) "parmread kontrol in cryst_theta"
795 !elwrite(iout,*) "parmread kontrol in cryst_theta"
797 athet(1,i,1,-1)=athet(1,i,1,1)
798 athet(2,i,1,-1)=athet(2,i,1,1)
799 bthet(1,i,1,-1)=-bthet(1,i,1,1)
800 bthet(2,i,1,-1)=-bthet(2,i,1,1)
801 athet(1,i,-1,1)=-athet(1,i,1,1)
802 athet(2,i,-1,1)=-athet(2,i,1,1)
803 bthet(1,i,-1,1)=bthet(1,i,1,1)
804 bthet(2,i,-1,1)=bthet(2,i,1,1)
806 !elwrite(iout,*) "parmread kontrol in cryst_theta"
809 athet(1,i,-1,-1)=athet(1,-i,1,1)
810 athet(2,i,-1,-1)=-athet(2,-i,1,1)
811 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
812 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
813 athet(1,i,-1,1)=athet(1,-i,1,1)
814 athet(2,i,-1,1)=-athet(2,-i,1,1)
815 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
816 bthet(2,i,-1,1)=bthet(2,-i,1,1)
817 athet(1,i,1,-1)=-athet(1,-i,1,1)
818 athet(2,i,1,-1)=athet(2,-i,1,1)
819 bthet(1,i,1,-1)=bthet(1,-i,1,1)
820 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
825 polthet(j,i)=polthet(j,-i)
828 gthet(j,i)=gthet(j,-i)
831 !elwrite(iout,*) "parmread kontrol in cryst_theta"
833 !elwrite(iout,*) "parmread kontrol in cryst_theta"
836 ! & 'Parameters of the virtual-bond valence angles:'
837 ! write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
838 ! & ' ATHETA0 ',' A1 ',' A2 ',
841 ! write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
842 ! & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
844 ! write (iout,'(/a/9x,5a/79(1h-))')
845 ! & 'Parameters of the expression for sigma(theta_c):',
846 ! & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
847 ! & ' ALPH3 ',' SIGMA0C '
849 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
850 ! & (polthet(j,i),j=0,3),sigc0(i)
852 ! write (iout,'(/a/9x,5a/79(1h-))')
853 ! & 'Parameters of the second gaussian:',
854 ! & ' THETA0 ',' SIGMA0 ',' G1 ',
857 ! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
858 ! & sig0(i),(gthet(j,i),j=1,3)
861 'Parameters of the virtual-bond valence angles:'
862 write (iout,'(/a/9x,5a/79(1h-))') &
863 'Coefficients of expansion',&
864 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
865 ' b1*10^1 ',' b2*10^1 '
867 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
868 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
869 (10*bthet(j,i,1,1),j=1,2)
871 write (iout,'(/a/9x,5a/79(1h-))') &
872 'Parameters of the expression for sigma(theta_c):',&
873 ' alpha0 ',' alph1 ',' alph2 ',&
874 ' alhp3 ',' sigma0c '
876 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
877 (polthet(j,i),j=0,3),sigc0(i)
879 write (iout,'(/a/9x,5a/79(1h-))') &
880 'Parameters of the second gaussian:',&
881 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
884 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
885 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
890 ! Read the parameters of Utheta determined from ab initio surfaces
891 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
893 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
894 ! write (iout,*) "tu dochodze"
895 IF (tor_mode.eq.0) THEN
896 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
897 read (ithep,*) nthetyp,ntheterm,ntheterm2,&
898 ntheterm3,nsingle,ndouble
899 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
901 !----------------------------------------------------
902 ! allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
903 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
904 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
905 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
906 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
907 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
908 !(maxtheterm,-maxthetyp1:maxthetyp1,&
909 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
910 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
911 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
912 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
913 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
914 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
915 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
916 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
917 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
918 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
919 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
920 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
921 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
922 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
923 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
924 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
925 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
928 read (ithep,*) (ithetyp(i),i=1,ntyp1)
930 ithetyp(i)=-ithetyp(-i)
932 ! write (iout,*) "tu dochodze"
933 aa0thet(:,:,:,:)=0.0d0
934 aathet(:,:,:,:,:)=0.0d0
935 bbthet(:,:,:,:,:,:)=0.0d0
936 ccthet(:,:,:,:,:,:)=0.0d0
937 ddthet(:,:,:,:,:,:)=0.0d0
938 eethet(:,:,:,:,:,:)=0.0d0
939 ffthet(:,:,:,:,:,:,:)=0.0d0
940 ggthet(:,:,:,:,:,:,:)=0.0d0
944 do j=-nthetyp,nthetyp
945 do k=-nthetyp,nthetyp
946 read (ithep,'(6a)') res1
947 read (ithep,*) aa0thet(i,j,k,iblock)
948 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
950 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
951 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
952 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
953 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
956 (((ffthet(llll,lll,ll,i,j,k,iblock),&
957 ffthet(lll,llll,ll,i,j,k,iblock),&
958 ggthet(llll,lll,ll,i,j,k,iblock),&
959 ggthet(lll,llll,ll,i,j,k,iblock),&
960 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
965 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
966 ! coefficients of theta-and-gamma-dependent terms are zero.
971 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
972 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
974 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
975 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
978 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
980 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
983 ! Substitution for D aminoacids from symmetry.
986 do j=-nthetyp,nthetyp
987 do k=-nthetyp,nthetyp
988 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
990 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
994 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
995 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
996 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
997 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1003 ffthet(llll,lll,ll,i,j,k,iblock)= &
1004 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1005 ffthet(lll,llll,ll,i,j,k,iblock)= &
1006 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1007 ggthet(llll,lll,ll,i,j,k,iblock)= &
1008 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1009 ggthet(lll,llll,ll,i,j,k,iblock)= &
1010 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1020 ! Control printout of the coefficients of virtual-bond-angle potentials
1024 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1028 write (iout,'(//4a)') &
1029 'Type ',onelett(i),onelett(j),onelett(k)
1030 write (iout,'(//a,10x,a)') " l","a[l]"
1031 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1032 write (iout,'(i2,1pe15.5)') &
1033 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1035 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
1036 "b",l,"c",l,"d",l,"e",l
1038 write (iout,'(i2,4(1pe15.5))') m,&
1039 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1040 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1044 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
1045 "f+",l,"f-",l,"g+",l,"g-",l
1048 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1049 ffthet(n,m,l,i,j,k,iblock),&
1050 ffthet(m,n,l,i,j,k,iblock),&
1051 ggthet(n,m,l,i,j,k,iblock),&
1052 ggthet(m,n,l,i,j,k,iblock)
1063 !C here will be the apropriate recalibrating for D-aminoacid
1064 read (ithep,*) nthetyp
1065 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1066 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1067 do i=-nthetyp+1,nthetyp-1
1068 read (ithep,*) nbend_kcc_Tb(i)
1069 do j=0,nbend_kcc_Tb(i)
1070 read (ithep,*) ijunk,v1bend_chyb(j,i)
1074 write (iout,'(a)') &
1075 "Parameters of the valence-only potentials"
1076 do i=-nthetyp+1,nthetyp-1
1077 write (iout,'(2a)') "Type ",toronelet(i)
1078 do k=0,nbend_kcc_Tb(i)
1079 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1086 !--------------- Reading theta parameters for nucleic acid-------
1087 read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,&
1088 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1089 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1090 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1091 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1092 nthetyp_nucl+1,nthetyp_nucl+1))
1093 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1094 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1095 nthetyp_nucl+1,nthetyp_nucl+1))
1096 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1097 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1098 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1099 nthetyp_nucl+1,nthetyp_nucl+1))
1100 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1101 nthetyp_nucl+1,nthetyp_nucl+1))
1102 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1103 nthetyp_nucl+1,nthetyp_nucl+1))
1104 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1105 nthetyp_nucl+1,nthetyp_nucl+1))
1106 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1107 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1108 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1109 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1110 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1111 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1113 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1114 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1116 read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1118 aa0thet_nucl(:,:,:)=0.0d0
1119 aathet_nucl(:,:,:,:)=0.0d0
1120 bbthet_nucl(:,:,:,:,:)=0.0d0
1121 ccthet_nucl(:,:,:,:,:)=0.0d0
1122 ddthet_nucl(:,:,:,:,:)=0.0d0
1123 eethet_nucl(:,:,:,:,:)=0.0d0
1124 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1125 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1130 read (ithep_nucl,'(3a)') t1,t2,t3
1131 read (ithep_nucl,*) aa0thet_nucl(i,j,k)
1132 read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1133 read (ithep_nucl,*) &
1134 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1135 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1136 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1137 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1138 read (ithep_nucl,*) &
1139 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1140 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1141 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1147 !-------------------------------------------
1148 allocate(nlob(ntyp1)) !(ntyp1)
1149 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1150 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1151 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1169 gaussc(l,k,j,i)=0.0D0
1177 ! Read the parameters of the probability distribution/energy expression
1178 ! of the side chains.
1181 !c write (iout,*) "tu dochodze",i
1182 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
1186 dsc_inv(i)=1.0D0/dsc(i)
1197 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
1198 censc(1,1,-i)=censc(1,1,i)
1199 censc(2,1,-i)=censc(2,1,i)
1200 censc(3,1,-i)=-censc(3,1,i)
1202 read (irotam,*) bsc(j,i)
1203 read (irotam,*) (censc(k,j,i),k=1,3),&
1204 ((blower(k,l,j),l=1,k),k=1,3)
1205 censc(1,j,-i)=censc(1,j,i)
1206 censc(2,j,-i)=censc(2,j,i)
1207 censc(3,j,-i)=-censc(3,j,i)
1208 ! BSC is amplitude of Gaussian
1215 akl=akl+blower(k,m,j)*blower(l,m,j)
1219 if (((k.eq.3).and.(l.ne.3)) &
1220 .or.((l.eq.3).and.(k.ne.3))) then
1221 gaussc(k,l,j,-i)=-akl
1222 gaussc(l,k,j,-i)=-akl
1224 gaussc(k,l,j,-i)=akl
1225 gaussc(l,k,j,-i)=akl
1234 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1237 if (nlobi.gt.0) then
1238 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1239 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1240 ! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1241 ! write (iout,'(a,f10.4,4(16x,f10.4))')
1242 ! & 'Center ',(bsc(j,i),j=1,nlobi)
1243 ! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
1244 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1245 'log h',(bsc(j,i),j=1,nlobi)
1246 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1247 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1248 ! write (iout,'(a)')
1254 ! blower(k,l,j)=gaussc(ind,j,i)
1259 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1260 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1267 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1268 ! added by Urszula Kozlowska 07/11/2007
1270 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1278 read(irotam,*) sc_parmin(j,i)
1284 !---------reading nucleic acid parameters for rotamers-------------------
1285 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1286 do i=1,ntyp_molec(2)
1287 read (irotam_nucl,*)
1289 read(irotam_nucl,*) sc_parmin_nucl(j,i)
1295 write (iout,*) "Base rotamer parameters"
1296 do i=1,ntyp_molec(2)
1297 write (iout,'(a)') restyp(i,2)
1298 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1303 read (ifourier,*) nloctyp
1304 !el write(iout,*)"nloctyp",nloctyp
1305 SPLIT_FOURIERTOR = nloctyp.lt.0
1306 nloctyp = iabs(nloctyp)
1308 if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1))
1309 print *,"shape",shape(itype2loc)
1310 allocate(iloctyp(-nloctyp:nloctyp))
1311 allocate(bnew1(3,2,-nloctyp:nloctyp))
1312 allocate(bnew2(3,2,-nloctyp:nloctyp))
1313 allocate(ccnew(3,2,-nloctyp:nloctyp))
1314 allocate(ddnew(3,2,-nloctyp:nloctyp))
1315 allocate(e0new(3,-nloctyp:nloctyp))
1316 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1317 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1318 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1319 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1320 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1321 allocate(e0newtor(3,-nloctyp:nloctyp))
1322 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1324 read (ifourier,*) (itype2loc(i),i=1,ntyp)
1325 read (ifourier,*) (iloctyp(i),i=0,nloctyp-1)
1326 itype2loc(ntyp1)=nloctyp
1327 iloctyp(nloctyp)=ntyp1
1329 itype2loc(-i)=-itype2loc(i)
1332 allocate(iloctyp(-nloctyp:nloctyp))
1333 allocate(itype2loc(-ntyp1:ntyp1))
1340 iloctyp(-i)=-iloctyp(i)
1342 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1343 !c write (iout,*) "nloctyp",nloctyp,
1344 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1345 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1346 !c write (iout,*) "nloctyp",nloctyp,
1347 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1350 !c write (iout,*) "NEWCORR",i
1354 read (ifourier,*) bnew1(ii,j,i)
1357 !c write (iout,*) "NEWCORR BNEW1"
1358 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1361 read (ifourier,*) bnew2(ii,j,i)
1364 !c write (iout,*) "NEWCORR BNEW2"
1365 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1367 read (ifourier,*) ccnew(kk,1,i)
1368 read (ifourier,*) ccnew(kk,2,i)
1370 !c write (iout,*) "NEWCORR CCNEW"
1371 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1373 read (ifourier,*) ddnew(kk,1,i)
1374 read (ifourier,*) ddnew(kk,2,i)
1376 !c write (iout,*) "NEWCORR DDNEW"
1377 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1381 read (ifourier,*) eenew(ii,jj,kk,i)
1385 !c write (iout,*) "NEWCORR EENEW1"
1386 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1388 read (ifourier,*) e0new(ii,i)
1390 !c write (iout,*) (e0new(ii,i),ii=1,3)
1392 !c write (iout,*) "NEWCORR EENEW"
1395 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1396 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1397 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1398 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1403 bnew1(ii,1,-i)= bnew1(ii,1,i)
1404 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1405 bnew2(ii,1,-i)= bnew2(ii,1,i)
1406 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1409 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1410 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1411 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1412 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1413 ccnew(ii,1,-i)=ccnew(ii,1,i)
1414 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1415 ddnew(ii,1,-i)=ddnew(ii,1,i)
1416 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1418 e0new(1,-i)= e0new(1,i)
1419 e0new(2,-i)=-e0new(2,i)
1420 e0new(3,-i)=-e0new(3,i)
1422 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1423 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1424 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1425 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1429 write (iout,'(a)') "Coefficients of the multibody terms"
1430 do i=-nloctyp+1,nloctyp-1
1431 write (iout,*) "Type: ",onelet(iloctyp(i))
1432 write (iout,*) "Coefficients of the expansion of B1"
1434 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1436 write (iout,*) "Coefficients of the expansion of B2"
1438 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1440 write (iout,*) "Coefficients of the expansion of C"
1441 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1442 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1443 write (iout,*) "Coefficients of the expansion of D"
1444 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1445 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1446 write (iout,*) "Coefficients of the expansion of E"
1447 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1450 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1455 IF (SPLIT_FOURIERTOR) THEN
1457 !c write (iout,*) "NEWCORR TOR",i
1461 read (ifourier,*) bnew1tor(ii,j,i)
1464 !c write (iout,*) "NEWCORR BNEW1 TOR"
1465 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1468 read (ifourier,*) bnew2tor(ii,j,i)
1471 !c write (iout,*) "NEWCORR BNEW2 TOR"
1472 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1474 read (ifourier,*) ccnewtor(kk,1,i)
1475 read (ifourier,*) ccnewtor(kk,2,i)
1477 !c write (iout,*) "NEWCORR CCNEW TOR"
1478 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1480 read (ifourier,*) ddnewtor(kk,1,i)
1481 read (ifourier,*) ddnewtor(kk,2,i)
1483 !c write (iout,*) "NEWCORR DDNEW TOR"
1484 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1488 read (ifourier,*) eenewtor(ii,jj,kk,i)
1492 !c write (iout,*) "NEWCORR EENEW1 TOR"
1493 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1495 read (ifourier,*) e0newtor(ii,i)
1497 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1499 !c write (iout,*) "NEWCORR EENEW TOR"
1502 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1503 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1504 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1505 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1510 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1511 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1512 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1513 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1516 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1517 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1518 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1519 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1521 e0newtor(1,-i)= e0newtor(1,i)
1522 e0newtor(2,-i)=-e0newtor(2,i)
1523 e0newtor(3,-i)=-e0newtor(3,i)
1525 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1526 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1527 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1528 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1532 write (iout,'(a)') &
1533 "Single-body coefficients of the torsional potentials"
1534 do i=-nloctyp+1,nloctyp-1
1535 write (iout,*) "Type: ",onelet(iloctyp(i))
1536 write (iout,*) "Coefficients of the expansion of B1tor"
1538 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1539 j,(bnew1tor(k,j,i),k=1,3)
1541 write (iout,*) "Coefficients of the expansion of B2tor"
1543 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1544 j,(bnew2tor(k,j,i),k=1,3)
1546 write (iout,*) "Coefficients of the expansion of Ctor"
1547 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1548 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1549 write (iout,*) "Coefficients of the expansion of Dtor"
1550 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1551 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1552 write (iout,*) "Coefficients of the expansion of Etor"
1553 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1556 write (iout,'(1hE,2i1,2f10.5)') &
1557 j,k,(eenewtor(l,j,k,i),l=1,2)
1563 do i=-nloctyp+1,nloctyp-1
1566 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1571 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1575 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1576 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1577 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1578 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1581 ENDIF !SPLIT_FOURIER_TOR
1583 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1584 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1585 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1586 allocate(b(13,-nloctyp-1:nloctyp+1))
1588 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1591 read (ifourier,*) (b(ii,i),ii=1,13)
1593 write (iout,*) 'Type ',onelet(iloctyp(i))
1594 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1602 CCold(1,1,i)= b(7,i)
1603 CCold(2,2,i)=-b(7,i)
1604 CCold(2,1,i)= b(9,i)
1605 CCold(1,2,i)= b(9,i)
1606 CCold(1,1,-i)= b(7,i)
1607 CCold(2,2,-i)=-b(7,i)
1608 CCold(2,1,-i)=-b(9,i)
1609 CCold(1,2,-i)=-b(9,i)
1610 DDold(1,1,i)= b(6,i)
1611 DDold(2,2,i)=-b(6,i)
1612 DDold(2,1,i)= b(8,i)
1613 DDold(1,2,i)= b(8,i)
1614 DDold(1,1,-i)= b(6,i)
1615 DDold(2,2,-i)=-b(6,i)
1616 DDold(2,1,-i)=-b(8,i)
1617 DDold(1,2,-i)=-b(8,i)
1618 EEold(1,1,i)= b(10,i)+b(11,i)
1619 EEold(2,2,i)=-b(10,i)+b(11,i)
1620 EEold(2,1,i)= b(12,i)-b(13,i)
1621 EEold(1,2,i)= b(12,i)+b(13,i)
1622 EEold(1,1,-i)= b(10,i)+b(11,i)
1623 EEold(2,2,-i)=-b(10,i)+b(11,i)
1624 EEold(2,1,-i)=-b(12,i)+b(13,i)
1625 EEold(1,2,-i)=-b(12,i)-b(13,i)
1626 write(iout,*) "TU DOCHODZE"
1632 "Coefficients of the cumulants (independent of valence angles)"
1633 do i=-nloctyp+1,nloctyp-1
1634 write (iout,*) 'Type ',onelet(iloctyp(i))
1636 write(iout,'(2f10.5)') B(3,i),B(5,i)
1638 write(iout,'(2f10.5)') B(2,i),B(4,i)
1641 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1645 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1649 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1656 ! Read torsional parameters in old format
1658 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1660 read (itorp,*) ntortyp,nterm_old
1661 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1662 read (itorp,*) (itortyp(i),i=1,ntyp)
1664 !el from energy module--------
1665 allocate(v1(nterm_old,ntortyp,ntortyp))
1666 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1667 !el---------------------------
1673 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1679 write (iout,'(/a/)') 'Torsional constants:'
1682 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1683 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1691 ! Read torsional parameters
1693 IF (TOR_MODE.eq.0) THEN
1695 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1697 read (itorp,*) ntortyp
1698 read (itorp,*) (itortyp(i),i=1,ntyp)
1699 write (iout,*) 'ntortyp',ntortyp
1701 !el from energy module---------
1702 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1703 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1705 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1706 allocate(vlor2(maxlor,ntortyp,ntortyp))
1707 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1708 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1710 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1711 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1712 !el---------------------------
1714 do i=-ntortyp,ntortyp
1715 do j=-ntortyp,ntortyp
1721 !el---------------------------
1725 itortyp(i)=-itortyp(-i)
1727 ! write (iout,*) 'ntortyp',ntortyp
1729 do j=-ntortyp+1,ntortyp-1
1730 read (itorp,*) nterm(i,j,iblock),&
1732 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1733 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1736 do k=1,nterm(i,j,iblock)
1737 read (itorp,*) kk,v1(k,i,j,iblock),&
1739 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1740 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1741 v0ij=v0ij+si*v1(k,i,j,iblock)
1744 do k=1,nlor(i,j,iblock)
1745 read (itorp,*) kk,vlor1(k,i,j),&
1746 vlor2(k,i,j),vlor3(k,i,j)
1747 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1750 v0(-i,-j,iblock)=v0ij
1757 write (iout,'(/a/)') 'Torsional constants:'
1760 write (iout,*) 'ityp',i,' jtyp',j
1761 write (iout,*) 'Fourier constants'
1762 do k=1,nterm(i,j,iblock)
1763 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1766 write (iout,*) 'Lorenz constants'
1767 do k=1,nlor(i,j,iblock)
1768 write (iout,'(3(1pe15.5))') &
1769 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1776 ! 6/23/01 Read parameters for double torsionals
1778 !el from energy module------------
1779 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1780 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1781 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1782 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1783 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1784 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1785 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1786 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1787 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1788 !---------------------------------
1792 do j=-ntortyp+1,ntortyp-1
1793 do k=-ntortyp+1,ntortyp-1
1794 read (itordp,'(3a1)') t1,t2,t3
1795 ! write (iout,*) "OK onelett",
1798 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1799 .or. t3.ne.toronelet(k)) then
1800 write (iout,*) "Error in double torsional parameter file",&
1803 call MPI_Finalize(Ierror)
1805 stop "Error in double torsional parameter file"
1807 read (itordp,*) ntermd_1(i,j,k,iblock),&
1808 ntermd_2(i,j,k,iblock)
1809 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1810 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1811 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1812 ntermd_1(i,j,k,iblock))
1813 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1814 ntermd_1(i,j,k,iblock))
1815 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1816 ntermd_1(i,j,k,iblock))
1817 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1818 ntermd_1(i,j,k,iblock))
1819 ! Martix of D parameters for one dimesional foureir series
1820 do l=1,ntermd_1(i,j,k,iblock)
1821 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1822 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1823 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1824 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1825 ! write(iout,*) "whcodze" ,
1826 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1828 read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1829 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1830 v2s(m,l,i,j,k,iblock),&
1831 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1832 ! Martix of D parameters for two dimesional fourier series
1833 do l=1,ntermd_2(i,j,k,iblock)
1835 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1836 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1837 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1838 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1847 write (iout,*) 'Constants for double torsionals'
1850 do j=-ntortyp+1,ntortyp-1
1851 do k=-ntortyp+1,ntortyp-1
1852 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1853 ' nsingle',ntermd_1(i,j,k,iblock),&
1854 ' ndouble',ntermd_2(i,j,k,iblock)
1856 write (iout,*) 'Single angles:'
1857 do l=1,ntermd_1(i,j,k,iblock)
1858 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1859 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1860 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1861 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1864 write (iout,*) 'Pairs of angles:'
1865 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1866 do l=1,ntermd_2(i,j,k,iblock)
1867 write (iout,'(i5,20f10.5)') &
1868 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1871 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1872 do l=1,ntermd_2(i,j,k,iblock)
1873 write (iout,'(i5,20f10.5)') &
1874 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1875 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1885 itype2loc(i)=itortyp(i)
1888 ELSE IF (TOR_MODE.eq.1) THEN
1890 !C read valence-torsional parameters
1891 read (itorp,*) ntortyp
1893 write (iout,*) "Valence-torsional parameters read in ntortyp",&
1895 read (itorp,*) (itortyp(i),i=1,ntyp)
1896 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1899 itype2loc(i)=itortyp(i)
1903 itortyp(i)=-itortyp(-i)
1905 do i=-ntortyp+1,ntortyp-1
1906 do j=-ntortyp+1,ntortyp-1
1907 !C first we read the cos and sin gamma parameters
1908 read (itorp,'(13x,a)') string
1909 write (iout,*) i,j,string
1911 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1912 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1913 do k=1,nterm_kcc(j,i)
1914 do l=1,nterm_kcc_Tb(j,i)
1915 do ll=1,nterm_kcc_Tb(j,i)
1916 read (itorp,*) ii,jj,kk, &
1917 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1925 !c AL 4/8/16: Calculate coefficients from one-body parameters
1927 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1928 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
1929 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
1930 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1931 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1934 print *,i,itortyp(i)
1935 itortyp(i)=itype2loc(i)
1938 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
1940 do i=-ntortyp+1,ntortyp-1
1941 do j=-ntortyp+1,ntortyp-1
1944 do k=1,nterm_kcc_Tb(j,i)
1945 do l=1,nterm_kcc_Tb(j,i)
1946 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
1947 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1948 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
1949 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1952 do k=1,nterm_kcc_Tb(j,i)
1953 do l=1,nterm_kcc_Tb(j,i)
1955 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1956 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1957 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1958 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1960 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1961 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1962 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1963 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1967 !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)
1971 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1975 if (tor_mode.gt.0 .and. lprint) then
1976 !c Print valence-torsional parameters
1977 write (iout,'(a)') &
1978 "Parameters of the valence-torsional potentials"
1979 do i=-ntortyp+1,ntortyp-1
1980 do j=-ntortyp+1,ntortyp-1
1981 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1982 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1983 do k=1,nterm_kcc(j,i)
1984 do l=1,nterm_kcc_Tb(j,i)
1985 do ll=1,nterm_kcc_Tb(j,i)
1986 write (iout,'(3i5,2f15.4)')&
1987 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1996 !elwrite(iout,*) "parmread kontrol sc-bb"
1997 ! Read of Side-chain backbone correlation parameters
1998 ! Modified 11 May 2012 by Adasko
2001 read (isccor,*) nsccortyp
2004 !c maxinter is maximum interaction sites
2005 !write(iout,*)"maxterm_sccor",maxterm_sccor
2006 !el from module energy-------------
2007 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2008 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2009 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2010 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2011 !-----------------------------------
2012 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2013 !-----------------------------------
2014 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2015 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2016 -nsccortyp:nsccortyp))
2017 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2018 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2019 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2020 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2021 !-----------------------------------
2022 do i=-nsccortyp,nsccortyp
2023 do j=-nsccortyp,nsccortyp
2027 !-----------------------------------
2029 read (isccor,*) (isccortyp(i),i=1,ntyp)
2031 isccortyp(i)=-isccortyp(-i)
2033 iscprol=isccortyp(20)
2034 ! write (iout,*) 'ntortyp',ntortyp
2036 !c maxinter is maximum interaction sites
2041 nterm_sccor(i,j),nlor_sccor(i,j)
2047 nterm_sccor(-i,j)=nterm_sccor(i,j)
2048 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2049 nterm_sccor(i,-j)=nterm_sccor(i,j)
2050 do k=1,nterm_sccor(i,j)
2051 read (isccor,*) kk,v1sccor(k,l,i,j),&
2053 if (j.eq.iscprol) then
2054 if (i.eq.isccortyp(10)) then
2055 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2056 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2058 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2059 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2060 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2061 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2062 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2063 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2064 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2065 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2068 if (i.eq.isccortyp(10)) then
2069 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2070 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2072 if (j.eq.isccortyp(10)) then
2073 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2074 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2076 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2077 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2078 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2079 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2080 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2081 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2085 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2086 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2087 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2088 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2091 do k=1,nlor_sccor(i,j)
2092 read (isccor,*) kk,vlor1sccor(k,i,j),&
2093 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2094 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2095 (1+vlor3sccor(k,i,j)**2)
2097 v0sccor(l,i,j)=v0ijsccor
2098 v0sccor(l,-i,j)=v0ijsccor1
2099 v0sccor(l,i,-j)=v0ijsccor2
2100 v0sccor(l,-i,-j)=v0ijsccor3
2106 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
2109 write (iout,*) 'ityp',i,' jtyp',j
2110 write (iout,*) 'Fourier constants'
2111 do k=1,nterm_sccor(i,j)
2112 write (iout,'(2(1pe15.5))') &
2113 (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
2115 write (iout,*) 'Lorenz constants'
2116 do k=1,nlor_sccor(i,j)
2117 write (iout,'(3(1pe15.5))') &
2118 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2125 ! Read electrostatic-interaction parameters
2128 write (iout,'(/a)') 'Electrostatic interaction constants:'
2129 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2130 'IT','JT','APP','BPP','AEL6','AEL3'
2132 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
2133 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
2134 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
2135 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
2140 app (i,j)=epp(i,j)*rri*rri
2141 bpp (i,j)=-2.0D0*epp(i,j)*rri
2142 ael6(i,j)=elpp6(i,j)*4.2D0**6
2143 ael3(i,j)=elpp3(i,j)*4.2D0**3
2144 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2149 ! Read side-chain interaction parameters.
2151 !el from module energy - COMMON.INTERACT-------
2152 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2153 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2154 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2155 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2156 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2157 allocate(epslip(ntyp,ntyp))
2168 !--------------------------------
2170 read (isidep,*) ipot,expon
2171 !el if (ipot.lt.1 .or. ipot.gt.5) then
2172 ! write (iout,'(2a)') 'Error while reading SC interaction',&
2173 ! 'potential file - unknown potential type.'
2177 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2178 ', exponents are ',expon,2*expon
2179 ! goto (10,20,30,30,40) ipot
2181 !----------------------- LJ potential ---------------------------------
2183 ! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2184 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2186 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2187 write (iout,'(a/)') 'The epsilon array:'
2188 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2189 write (iout,'(/a)') 'One-body parameters:'
2190 write (iout,'(a,4x,a)') 'residue','sigma'
2191 write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
2194 !----------------------- LJK potential --------------------------------
2196 ! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2197 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2198 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2200 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2201 write (iout,'(a/)') 'The epsilon array:'
2202 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2203 write (iout,'(/a)') 'One-body parameters:'
2204 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2205 write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2209 !---------------------- GB or BP potential -----------------------------
2212 if (scelemode.eq.0) then
2214 read (isidep,*)(eps(i,j),j=i,ntyp)
2216 read (isidep,*)(sigma0(i),i=1,ntyp)
2217 read (isidep,*)(sigii(i),i=1,ntyp)
2218 read (isidep,*)(chip(i),i=1,ntyp)
2219 read (isidep,*)(alp(i),i=1,ntyp)
2221 read (isidep,*)(epslip(i,j),j=i,ntyp)
2223 ! For the GB potential convert sigma'**2 into chi'
2226 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2230 write (iout,'(/a/)') 'Parameters of the BP potential:'
2231 write (iout,'(a/)') 'The epsilon array:'
2232 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2233 write (iout,'(/a)') 'One-body parameters:'
2234 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2236 write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
2237 chip(i),alp(i),i=1,ntyp)
2240 allocate(icharge(ntyp1))
2241 ! print *,ntyp,icharge(i)
2243 read (isidep,*) (icharge(i),i=1,ntyp)
2244 print *,ntyp,icharge(i)
2245 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2246 write (2,*) "icharge",(icharge(i),i=1,ntyp)
2247 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2248 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2249 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2250 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2251 allocate(epsintab(ntyp,ntyp))
2252 allocate(dtail(2,ntyp,ntyp))
2253 print *,"control line 1"
2254 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2255 allocate(wqdip(2,ntyp,ntyp))
2256 allocate(wstate(4,ntyp,ntyp))
2257 allocate(dhead(2,2,ntyp,ntyp))
2258 allocate(nstate(ntyp,ntyp))
2259 allocate(debaykap(ntyp,ntyp))
2260 print *,"control line 2"
2261 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2262 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2266 ! write (*,*) "Im in ALAB", i, " ", j
2268 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2269 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2270 chis(i,j),chis(j,i), &
2271 nstate(i,j),(wstate(k,i,j),k=1,4), &
2272 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2273 dtail(1,i,j),dtail(2,i,j), &
2274 epshead(i,j),sig0head(i,j), &
2275 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2276 alphapol(i,j),alphapol(j,i), &
2277 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
2278 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2284 sigma(i,j) = sigma(j,i)
2285 sigmap1(i,j)=sigmap1(j,i)
2286 sigmap2(i,j)=sigmap2(j,i)
2287 sigiso1(i,j)=sigiso1(j,i)
2288 sigiso2(i,j)=sigiso2(j,i)
2289 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2290 nstate(i,j) = nstate(j,i)
2291 dtail(1,i,j) = dtail(1,j,i)
2292 dtail(2,i,j) = dtail(2,j,i)
2294 alphasur(k,i,j) = alphasur(k,j,i)
2295 wstate(k,i,j) = wstate(k,j,i)
2296 alphiso(k,i,j) = alphiso(k,j,i)
2299 dhead(2,1,i,j) = dhead(1,1,j,i)
2300 dhead(2,2,i,j) = dhead(1,2,j,i)
2301 dhead(1,1,i,j) = dhead(2,1,j,i)
2302 dhead(1,2,i,j) = dhead(2,2,j,i)
2304 epshead(i,j) = epshead(j,i)
2305 sig0head(i,j) = sig0head(j,i)
2308 wqdip(k,i,j) = wqdip(k,j,i)
2311 wquad(i,j) = wquad(j,i)
2312 epsintab(i,j) = epsintab(j,i)
2313 debaykap(i,j)=debaykap(j,i)
2314 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2320 !--------------------- GBV potential -----------------------------------
2322 ! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2323 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2324 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2325 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2327 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2328 write (iout,'(a/)') 'The epsilon array:'
2329 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2330 write (iout,'(/a)') 'One-body parameters:'
2331 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2332 's||/s_|_^2',' chip ',' alph '
2333 write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2334 sigii(i),chip(i),alp(i),i=1,ntyp)
2337 write (iout,'(2a)') 'Error while reading SC interaction',&
2338 'potential file - unknown potential type.'
2344 !-----------------------------------------------------------------------
2345 ! Calculate the "working" parameters of SC interactions.
2347 !el from module energy - COMMON.INTERACT-------
2348 ! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2349 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1))
2350 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) !(ntyp,ntyp)
2351 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2352 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2353 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2354 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2356 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2364 if (scelemode.eq.0) then
2371 !--------------------------------
2376 epslip(i,j)=epslip(j,i)
2379 if (scelemode.eq.0) then
2382 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2383 sigma(j,i)=sigma(i,j)
2384 rs0(i,j)=dwa16*sigma(i,j)
2389 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2390 'Working parameters of the SC interactions:',&
2391 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2396 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2406 sigeps=dsign(1.0D0,epsij)
2408 aa_aq(i,j)=epsij*rrij*rrij
2409 bb_aq(i,j)=-sigeps*epsij*rrij
2410 aa_aq(j,i)=aa_aq(i,j)
2411 bb_aq(j,i)=bb_aq(i,j)
2412 epsijlip=epslip(i,j)
2413 sigeps=dsign(1.0D0,epsijlip)
2414 epsijlip=dabs(epsijlip)
2415 aa_lip(i,j)=epsijlip*rrij*rrij
2416 bb_lip(i,j)=-sigeps*epsijlip*rrij
2417 aa_lip(j,i)=aa_lip(i,j)
2418 bb_lip(j,i)=bb_lip(i,j)
2419 if ((ipot.gt.2).and. (scelemode.eq.0))then
2420 sigt1sq=sigma0(i)**2
2421 sigt2sq=sigma0(j)**2
2424 ratsig1=sigt2sq/sigt1sq
2425 ratsig2=1.0D0/ratsig1
2426 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2427 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2428 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2432 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2433 sigmaii(i,j)=rsum_max
2434 sigmaii(j,i)=rsum_max
2436 ! sigmaii(i,j)=r0(i,j)
2437 ! sigmaii(j,i)=r0(i,j)
2439 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2440 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2441 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2442 augm(i,j)=epsij*r_augm**(2*expon)
2443 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2450 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2451 restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2452 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2456 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2457 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2458 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2459 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2460 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2461 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2462 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2463 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2464 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2465 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2466 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2467 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2468 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2469 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2470 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2471 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2480 read (isidep_nucl,*) ipot_nucl
2481 ! print *,"TU?!",ipot_nucl
2482 if (ipot_nucl.eq.1) then
2483 do i=1,ntyp_molec(2)
2484 do j=i,ntyp_molec(2)
2485 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2486 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2490 do i=1,ntyp_molec(2)
2491 do j=i,ntyp_molec(2)
2492 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2493 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2494 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2498 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2499 do i=1,ntyp_molec(2)
2500 do j=i,ntyp_molec(2)
2501 rrij=sigma_nucl(i,j)
2505 epsij=4*eps_nucl(i,j)
2506 sigeps=dsign(1.0D0,epsij)
2508 aa_nucl(i,j)=epsij*rrij*rrij
2509 bb_nucl(i,j)=-sigeps*epsij*rrij
2510 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2511 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2512 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2513 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2514 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2515 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2518 aa_nucl(i,j)=aa_nucl(j,i)
2519 bb_nucl(i,j)=bb_nucl(j,i)
2520 ael3_nucl(i,j)=ael3_nucl(j,i)
2521 ael6_nucl(i,j)=ael6_nucl(j,i)
2522 ael63_nucl(i,j)=ael63_nucl(j,i)
2523 ael32_nucl(i,j)=ael32_nucl(j,i)
2524 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2525 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2526 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2527 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2528 eps_nucl(i,j)=eps_nucl(j,i)
2529 sigma_nucl(i,j)=sigma_nucl(j,i)
2530 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2534 write(iout,*) "tube param"
2535 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2536 ccavtubpep,dcavtubpep,tubetranenepep
2537 sigmapeptube=sigmapeptube**6
2538 sigeps=dsign(1.0D0,epspeptube)
2539 epspeptube=dabs(epspeptube)
2540 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2541 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2542 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2544 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2545 ccavtub(i),dcavtub(i),tubetranene(i)
2546 sigmasctube=sigmasctube**6
2547 sigeps=dsign(1.0D0,epssctube)
2548 epssctube=dabs(epssctube)
2549 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2550 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2551 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2553 !-----------------READING SC BASE POTENTIALS-----------------------------
2554 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2555 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2556 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2557 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2558 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2559 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2560 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2561 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2562 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2563 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2564 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2565 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2566 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2567 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2568 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2569 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2572 do i=1,ntyp_molec(1)
2573 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2574 write (*,*) "Im in ", i, " ", j
2575 read(isidep_scbase,*) &
2576 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2577 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2578 write(*,*) "eps",eps_scbase(i,j)
2579 read(isidep_scbase,*) &
2580 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2581 chis_scbase(i,j,1),chis_scbase(i,j,2)
2582 read(isidep_scbase,*) &
2583 dhead_scbasei(i,j), &
2584 dhead_scbasej(i,j), &
2585 rborn_scbasei(i,j),rborn_scbasej(i,j)
2586 read(isidep_scbase,*) &
2587 (wdipdip_scbase(k,i,j),k=1,3), &
2588 (wqdip_scbase(k,i,j),k=1,2)
2589 read(isidep_scbase,*) &
2590 alphapol_scbase(i,j), &
2591 epsintab_scbase(i,j)
2594 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2595 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2597 do i=1,ntyp_molec(1)
2598 do j=1,ntyp_molec(2)-1
2599 epsij=eps_scbase(i,j)
2600 rrij=sigma_scbase(i,j)
2605 sigeps=dsign(1.0D0,epsij)
2607 aa_scbase(i,j)=epsij*rrij*rrij
2608 bb_scbase(i,j)=-sigeps*epsij*rrij
2611 !-----------------READING PEP BASE POTENTIALS-------------------
2612 allocate(eps_pepbase(ntyp_molec(2)))
2613 allocate(sigma_pepbase(ntyp_molec(2)))
2614 allocate(chi_pepbase(ntyp_molec(2),2))
2615 allocate(chipp_pepbase(ntyp_molec(2),2))
2616 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2617 allocate(sigmap1_pepbase(ntyp_molec(2)))
2618 allocate(sigmap2_pepbase(ntyp_molec(2)))
2619 allocate(chis_pepbase(ntyp_molec(2),2))
2620 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2623 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2624 write (*,*) "Im in ", i, " ", j
2625 read(isidep_pepbase,*) &
2626 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2627 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2628 write(*,*) "eps",eps_pepbase(j)
2629 read(isidep_pepbase,*) &
2630 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2631 chis_pepbase(j,1),chis_pepbase(j,2)
2632 read(isidep_pepbase,*) &
2633 (wdipdip_pepbase(k,j),k=1,3)
2635 allocate(aa_pepbase(ntyp_molec(2)))
2636 allocate(bb_pepbase(ntyp_molec(2)))
2638 do j=1,ntyp_molec(2)-1
2639 epsij=eps_pepbase(j)
2640 rrij=sigma_pepbase(j)
2645 sigeps=dsign(1.0D0,epsij)
2647 aa_pepbase(j)=epsij*rrij*rrij
2648 bb_pepbase(j)=-sigeps*epsij*rrij
2650 !--------------READING SC PHOSPHATE-------------------------------------
2651 !--------------READING SC PHOSPHATE-------------------------------------
2652 allocate(eps_scpho(ntyp_molec(1)))
2653 allocate(sigma_scpho(ntyp_molec(1)))
2654 allocate(chi_scpho(ntyp_molec(1),2))
2655 allocate(chipp_scpho(ntyp_molec(1),2))
2656 allocate(alphasur_scpho(4,ntyp_molec(1)))
2657 allocate(sigmap1_scpho(ntyp_molec(1)))
2658 allocate(sigmap2_scpho(ntyp_molec(1)))
2659 allocate(chis_scpho(ntyp_molec(1),2))
2660 allocate(wqq_scpho(ntyp_molec(1)))
2661 allocate(wqdip_scpho(2,ntyp_molec(1)))
2662 allocate(alphapol_scpho(ntyp_molec(1)))
2663 allocate(epsintab_scpho(ntyp_molec(1)))
2664 allocate(dhead_scphoi(ntyp_molec(1)))
2665 allocate(rborn_scphoi(ntyp_molec(1)))
2666 allocate(rborn_scphoj(ntyp_molec(1)))
2667 allocate(alphi_scpho(ntyp_molec(1)))
2671 do j=1,ntyp_molec(1) ! without U then we will take T for U
2672 write (*,*) "Im in scpho ", i, " ", j
2673 read(isidep_scpho,*) &
2674 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2675 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2676 write(*,*) "eps",eps_scpho(j)
2677 read(isidep_scpho,*) &
2678 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2679 chis_scpho(j,1),chis_scpho(j,2)
2680 read(isidep_scpho,*) &
2681 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2682 read(isidep_scpho,*) &
2683 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2687 allocate(aa_scpho(ntyp_molec(1)))
2688 allocate(bb_scpho(ntyp_molec(1)))
2689 do j=1,ntyp_molec(1)
2696 sigeps=dsign(1.0D0,epsij)
2698 aa_scpho(j)=epsij*rrij*rrij
2699 bb_scpho(j)=-sigeps*epsij*rrij
2703 read(isidep_peppho,*) &
2704 eps_peppho,sigma_peppho
2705 read(isidep_peppho,*) &
2706 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2707 read(isidep_peppho,*) &
2708 (wqdip_peppho(k),k=1,2)
2716 sigeps=dsign(1.0D0,epsij)
2718 aa_peppho=epsij*rrij*rrij
2719 bb_peppho=-sigeps*epsij*rrij
2722 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2730 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
2731 allocate(alphapolcat2(ntyp,ntyp))
2732 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
2733 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
2734 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
2735 allocate(epsintabcat(ntyp,ntyp))
2736 allocate(dtailcat(2,ntyp,ntyp))
2737 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
2738 allocate(wqdipcat(2,ntyp,ntyp))
2739 allocate(wstatecat(4,ntyp,ntyp))
2740 allocate(dheadcat(2,2,ntyp,ntyp))
2741 allocate(nstatecat(ntyp,ntyp))
2742 allocate(debaykapcat(ntyp,ntyp))
2744 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
2745 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
2746 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
2747 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
2748 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
2751 allocate (ichargecat(ntyp_molec(5)))
2752 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
2753 if (oldion.eq.0) then
2754 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
2755 allocate(icharge(1:ntyp1))
2756 read(iion,*) (icharge(i),i=1,ntyp)
2761 do i=1,ntyp_molec(5)
2762 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
2763 print *,msc(i,5),restok(i,5)
2767 do j=1,ntyp_molec(5)
2769 ! do j=1,ntyp_molec(5)
2770 ! write (*,*) "Im in ALAB", i, " ", j
2772 epscat(i,j),sigmacat(i,j), &
2773 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
2774 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
2776 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
2777 ! chiscat(i,j),chiscat(j,i), &
2778 chis1cat(i,j),chis2cat(i,j), &
2780 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2781 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
2782 dtailcat(1,i,j),dtailcat(2,i,j), &
2783 epsheadcat(i,j),sig0headcat(i,j), &
2785 ! rborncat(i,j),rborncat(j,i),&
2786 rborn1cat(i,j),rborn2cat(i,j),&
2787 (wqdipcat(k,i,j),k=1,2), &
2788 alphapolcat(i,j),alphapolcat2(j,i), &
2789 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
2790 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2793 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
2795 do j=1,ntyp_molec(5)
2799 sigeps=dsign(1.0D0,epsij)
2801 aa_aq_cat(i,j)=epsij*rrij*rrij
2802 bb_aq_cat(i,j)=-sigeps*epsij*rrij
2806 do j=1,ntyp_molec(5)
2808 write (iout,*) 'i= ', i, ' j= ', j
2809 write (iout,*) 'epsi0= ', epscat(i,j)
2810 write (iout,*) 'sigma0= ', sigmacat(i,j)
2811 write (iout,*) 'chi1= ', chi1cat(i,j)
2812 write (iout,*) 'chi1= ', chi2cat(i,j)
2813 write (iout,*) 'chip1= ', chipp1cat(1,j)
2814 write (iout,*) 'chip2= ', chipp2cat(1,j)
2815 write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
2816 write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
2817 write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
2818 write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
2819 write (iout,*) 'sig1= ', sigmap1cat(1,j)
2820 write (iout,*) 'sig2= ', sigmap2cat(1,j)
2821 write (iout,*) 'chis1= ', chis1cat(1,j)
2822 write (iout,*) 'chis1= ', chis2cat(1,j)
2823 write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
2824 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
2825 write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
2826 write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
2827 write (iout,*) 'a1= ', rborn1cat(i,j)
2828 write (iout,*) 'a2= ', rborn2cat(i,j)
2829 write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
2830 write (iout,*) 'alphapol1= ', alphapolcat(1,j)
2831 write (iout,*) 'alphapol2= ', alphapolcat(j,1)
2832 write (iout,*) 'w1= ', wqdipcat(1,i,j)
2833 write (iout,*) 'w2= ', wqdipcat(2,i,j)
2834 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j)
2837 If ((i.eq.1).and.(j.eq.27)) then
2838 write (iout,*) 'SEP'
2839 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
2840 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
2850 ! Define the SC-p interaction constants
2859 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2860 read (itorp_nucl,*) ntortyp_nucl
2861 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2862 !el from energy module---------
2863 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2864 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2866 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2867 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2868 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2869 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2871 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2872 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2873 !el---------------------------
2876 !el--------------------
2877 read (itorp_nucl,*) &
2878 (itortyp_nucl(i),i=1,ntyp_molec(2))
2879 ! print *,itortyp_nucl(:)
2880 !c write (iout,*) 'ntortyp',ntortyp
2883 read (itorp_nucl,*) nterm_nucl(i,j),nlor_nucl(i,j)
2884 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2887 do k=1,nterm_nucl(i,j)
2888 read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2889 v0ij=v0ij+si*v1_nucl(k,i,j)
2892 do k=1,nlor_nucl(i,j)
2893 read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),&
2894 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2895 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2902 !elwrite(iout,*) "parmread kontrol before oldscp"
2904 ! Define the SC-p interaction constants
2908 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2910 ! aad(i,1)=0.3D0*4.0D0**12
2911 ! Following line for constants currently implemented
2912 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2913 aad(i,1)=1.5D0*4.0D0**12
2914 ! aad(i,1)=0.17D0*5.6D0**12
2916 ! "Soft" SC-p repulsion
2918 ! Following line for constants currently implemented
2919 ! aad(i,1)=0.3D0*4.0D0**6
2920 ! "Hard" SC-p repulsion
2921 bad(i,1)=3.0D0*4.0D0**6
2922 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2931 ! 8/9/01 Read the SC-p interaction constants from file
2934 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
2937 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2938 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2939 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2940 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2944 write (iout,*) "Parameters of SC-p interactions:"
2946 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2947 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2951 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2953 do i=1,ntyp_molec(2)
2954 read (iscpp_nucl,*) eps_scp_nucl(i),rscp_nucl(i)
2956 do i=1,ntyp_molec(2)
2957 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2958 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2960 r0pp=1.12246204830937298142*5.16158
2966 ! Define the constants of the disulfide bridge
2970 ! Old arbitrary potential - commented out.
2975 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2976 ! energy surface of diethyl disulfide.
2977 ! A. Liwo and U. Kozlowska, 11/24/03
2989 ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
2990 Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
2991 akcm=akcm/wsc*ssscale
2992 akth=akth/wsc*ssscale
2993 akct=akct/wsc*ssscale
2994 v1ss=v1ss/wsc*ssscale
2995 v2ss=v2ss/wsc*ssscale
2996 v3ss=v3ss/wsc*ssscale
2998 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
3002 write (iout,'(/a)') "Disulfide bridge parameters:"
3003 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3004 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3005 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3006 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3010 end subroutine parmread
3012 !-----------------------------------------------------------------------------
3014 !-----------------------------------------------------------------------------
3015 subroutine mygetenv(string,var)
3019 ! This subroutine passes the environmental variables to FORTRAN program.
3020 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
3021 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
3022 ! reads the environmental variables from $HOME/.env
3024 ! Usage: As for the standard FORTRAN GETENV subroutine.
3026 ! Purpose: some versions/installations of MPI do not transfer the environmental
3027 ! variables to slave processors, if these variables are set in the shell script
3028 ! from which mpirun is called.
3037 character*(*) :: string,var
3038 #if defined(MYGETENV) && defined(MPI)
3039 ! include "DIMENSIONS.ZSCOPT"
3041 ! include "COMMON.MPI"
3042 !el character*360 ucase
3044 character(len=360) :: string1(360),karta
3045 character(len=240) :: home
3048 call getenv("HOME",home)
3049 open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
3051 read (99,end=111,err=111,'(a)') karta
3055 call split_string(karta,string1,80,n)
3056 if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
3057 string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
3059 print *,"Processor",me,": ",var(:ilen(var)),&
3060 " assigned to ",string(:ilen(string))
3065 111 print *,"Environment variable ",string(:ilen(string))," not set."
3068 112 print *,"Error opening environment file!"
3070 call getenv(string,var)
3073 end subroutine mygetenv
3074 !-----------------------------------------------------------------------------
3076 !-----------------------------------------------------------------------------
3077 subroutine read_general_data(*)
3079 use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
3080 scelemode,TUBEmode,tor_mode,energy_dec
3082 use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss
3083 use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
3084 bordtubebot,tubebufthick,buftubebot,buftubetop
3086 ! include "DIMENSIONS"
3087 ! include "DIMENSIONS.ZSCOPT"
3088 ! include "DIMENSIONS.FREE"
3089 ! include "COMMON.TORSION"
3090 ! include "COMMON.INTERACT"
3091 ! include "COMMON.IOUNITS"
3092 ! include "COMMON.TIME1"
3093 ! include "COMMON.PROT"
3094 ! include "COMMON.PROTFILES"
3095 ! include "COMMON.CHAIN"
3096 ! include "COMMON.NAMES"
3097 ! include "COMMON.FFIELD"
3098 ! include "COMMON.ENEPS"
3099 ! include "COMMON.WEIGHTS"
3100 ! include "COMMON.FREE"
3101 ! include "COMMON.CONTROL"
3102 ! include "COMMON.ENERGIES"
3103 character(len=800) :: controlcard
3104 integer :: i,j,k,ii,n_ene_found
3105 integer :: ind,itype1,itype2,itypf,itypsc,itypp
3108 !el character*16 ucase
3109 character(len=16) :: key
3111 call card_concat(controlcard,.true.)
3112 call readi(controlcard,"N_ENE",n_eneW,max_eneW)
3113 if (n_eneW.gt.max_eneW) then
3114 write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
3118 call readi(controlcard,"NPARMSET",nparmset,1)
3119 !elwrite(iout,*)"in read_gen data"
3120 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
3121 call readi(controlcard,"IPARMPRINT",iparmprint,1)
3122 write (iout,*) "PARMPRINT",iparmprint
3123 if (nparmset.gt.max_parm) then
3124 write (iout,*) "Error: parameter out of range: NPARMSET",&
3128 !elwrite(iout,*)"in read_gen data"
3129 call readi(controlcard,"MAXIT",maxit,5000)
3130 call reada(controlcard,"FIMIN",fimin,1.0d-3)
3131 call readi(controlcard,"ENSEMBLES",ensembles,0)
3132 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
3133 write (iout,*) "Number of energy parameter sets",nparmset
3134 allocate(isampl(nparmset))
3135 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
3136 write (iout,*) "MaxSlice",MaxSlice
3137 call readi(controlcard,"NSLICE",nslice,1)
3138 !elwrite(iout,*)"in read_gen data"
3140 if (nslice.gt.MaxSlice) then
3141 write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
3145 write (iout,*) "Frequency of storing conformations",&
3146 (isampl(i),i=1,nparmset)
3147 write (iout,*) "Maxit",maxit," Fimin",fimin
3148 call readi(controlcard,"NQ",nQ,1)
3149 if (nQ.gt.MaxQ) then
3150 write (iout,*) "Error: parameter out of range: NQ",nq,&
3155 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
3156 call reada(controlcard,"DELTA",delta,1.0d-2)
3157 call readi(controlcard,"EINICHECK",einicheck,2)
3158 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
3159 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
3160 call readi(controlcard,"RESCALE",rescale_modeW,1)
3161 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3162 call reada(controlcard,'BOXY',boxysize,100.0d0)
3163 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3164 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3165 call readi(controlcard,"SCELEMODE",scelemode,0)
3166 call readi(controlcard,"OLDION",oldion,0)
3167 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
3168 print *,"SCELE",scelemode
3169 call readi(controlcard,'TORMODE',tor_mode,0)
3170 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3171 write(iout,*) "torsional and valence angle mode",tor_mode
3173 call readi(controlcard,'TUBEMOD',tubemode,0)
3175 if (TUBEmode.gt.0) then
3176 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3177 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3178 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3179 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3180 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3181 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3182 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3183 buftubebot=bordtubebot+tubebufthick
3184 buftubetop=bordtubetop-tubebufthick
3186 ions=index(controlcard,"IONS").gt.0
3187 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
3188 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3189 write(iout,*) "R_CUT_ELE=",r_cut_ele
3190 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
3191 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
3192 call readi(controlcard,'SYM',symetr,1)
3193 write (iout,*) "DISTCHAINMAX",distchainmax
3194 write (iout,*) "delta",delta
3195 write (iout,*) "einicheck",einicheck
3196 write (iout,*) "rescale_mode",rescale_modeW
3198 bxfile=index(controlcard,"BXFILE").gt.0
3199 cxfile=index(controlcard,"CXFILE").gt.0
3200 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
3202 histfile=index(controlcard,"HISTFILE").gt.0
3203 histout=index(controlcard,"HISTOUT").gt.0
3204 entfile=index(controlcard,"ENTFILE").gt.0
3205 zscfile=index(controlcard,"ZSCFILE").gt.0
3206 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
3207 write (iout,*) "with_dihed_constr ",with_dihed_constr
3208 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3210 end subroutine read_general_data
3211 !------------------------------------------------------------------------------
3212 subroutine read_efree(*)
3214 ! Read molecular data
3217 ! include 'DIMENSIONS'
3218 ! include 'DIMENSIONS.ZSCOPT'
3219 ! include 'DIMENSIONS.COMPAR'
3220 ! include 'DIMENSIONS.FREE'
3221 ! include 'COMMON.IOUNITS'
3222 ! include 'COMMON.TIME1'
3223 ! include 'COMMON.SBRIDGE'
3224 ! include 'COMMON.CONTROL'
3225 ! include 'COMMON.CHAIN'
3226 ! include 'COMMON.HEADER'
3227 ! include 'COMMON.GEO'
3228 ! include 'COMMON.FREE'
3229 character(len=320) :: controlcard !,ucase
3230 integer :: iparm,ib,i,j,npars
3240 ! call alloc_wham_arrays
3241 ! allocate(nT_h(nParmSet))
3242 ! allocate(replica(nParmSet))
3243 ! allocate(umbrella(nParmSet))
3244 ! allocate(read_iset(nParmSet))
3245 ! allocate(nT_h(nParmSet))
3249 call card_concat(controlcard,.true.)
3250 call readi(controlcard,'NT',nT_h(iparm),1)
3251 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
3253 if (nT_h(iparm).gt.MaxT_h) then
3254 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
3258 replica(iparm)=index(controlcard,"REPLICA").gt.0
3259 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
3260 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
3261 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
3262 replica(iparm)," umbrella ",umbrella(iparm),&
3263 " read_iset",read_iset(iparm)
3266 call card_concat(controlcard,.true.)
3267 call readi(controlcard,'NR',nR(ib,iparm),1)
3268 if (umbrella(iparm)) then
3271 nRR(ib,iparm)=nR(ib,iparm)
3273 if (nR(ib,iparm).gt.MaxR) then
3274 write (iout,*) "Error: parameter out of range: NR",&
3278 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
3279 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
3280 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
3283 call card_concat(controlcard,.true.)
3284 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
3286 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
3291 write (iout,*) "ib",ib," beta_h",&
3292 1.0d0/(0.001987*beta_h(ib,iparm))
3293 write (iout,*) "nR",nR(ib,iparm)
3294 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
3296 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
3297 "q0",(q0(j,i,ib,iparm),j=1,nQ)
3309 nR(ib,iparm)=nR(ib,1)
3310 if (umbrella(iparm)) then
3313 nRR(ib,iparm)=nR(ib,1)
3315 beta_h(ib,iparm)=beta_h(ib,1)
3317 f(i,ib,iparm)=f(i,ib,1)
3319 KH(j,i,ib,iparm)=KH(j,i,ib,1)
3320 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
3323 replica(iparm)=replica(1)
3324 umbrella(iparm)=umbrella(1)
3325 read_iset(iparm)=read_iset(1)
3332 end subroutine read_efree
3333 !-----------------------------------------------------------------------------
3334 subroutine read_protein_data(*)
3336 ! include "DIMENSIONS"
3337 ! include "DIMENSIONS.ZSCOPT"
3338 ! include "DIMENSIONS.FREE"
3342 integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
3343 ! include "COMMON.MPI"
3345 ! include "COMMON.CHAIN"
3346 ! include "COMMON.IOUNITS"
3347 ! include "COMMON.PROT"
3348 ! include "COMMON.PROTFILES"
3349 ! include "COMMON.NAMES"
3350 ! include "COMMON.FREE"
3351 ! include "COMMON.OBCINKA"
3352 character(len=64) :: nazwa
3353 character(len=16000) :: controlcard
3354 integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
3355 !el external ilen,iroof
3364 ! Read names of files with conformation data.
3365 if (replica(iparm)) then
3371 do ii=1,nRR(ib,iparm)
3372 write (iout,*) "Parameter set",iparm," temperature",ib,&
3375 call card_concat(controlcard,.true.)
3376 write (iout,*) controlcard(:ilen(controlcard))
3377 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
3378 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
3379 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
3380 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
3381 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
3382 maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
3383 call reada(controlcard,"TIME_START",&
3384 time_start_collect(ii,ib,iparm),0.0d0)
3385 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
3387 write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
3388 " rec_end",rec_end(ii,ib,iparm)
3389 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
3390 " time_end",time_end_collect(ii,ib,iparm)
3392 if (replica(iparm)) then
3393 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
3394 write (iout,*) "Number of trajectories",totraj(ii,iparm)
3397 if (nfile_bin(ii,ib,iparm).lt.2 &
3398 .and. nfile_asc(ii,ib,iparm).eq.0 &
3399 .and. nfile_cx(ii,ib,iparm).eq.0) then
3400 write (iout,*) "Error - no action specified!"
3403 if (nfile_bin(ii,ib,iparm).gt.0) then
3404 call card_concat(controlcard,.false.)
3405 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
3406 maxfile_prot,nfile_bin(ii,ib,iparm))
3408 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
3409 write(iout,*) (protfiles(i,1,ii,ib,iparm),&
3410 i=1,nfile_bin(ii,ib,iparm))
3413 if (nfile_asc(ii,ib,iparm).gt.0) then
3414 call card_concat(controlcard,.false.)
3415 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3416 maxfile_prot,nfile_asc(ii,ib,iparm))
3418 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
3419 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3420 i=1,nfile_asc(ii,ib,iparm))
3422 else if (nfile_cx(ii,ib,iparm).gt.0) then
3423 call card_concat(controlcard,.false.)
3424 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3425 maxfile_prot,nfile_cx(ii,ib,iparm))
3427 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
3428 write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3429 i=1,nfile_cx(ii,ib,iparm))
3438 end subroutine read_protein_data
3439 !-------------------------------------------------------------------------------
3440 subroutine readsss(rekord,lancuch,wartosc,default)
3442 character*(*) :: rekord,lancuch,wartosc,default
3443 character(len=80) :: aux
3444 integer :: lenlan,lenrec,iread,ireade
3448 lenlan=ilen(lancuch)
3450 iread=index(rekord,lancuch(:lenlan)//"=")
3451 ! print *,"rekord",rekord," lancuch",lancuch
3452 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3453 if (iread.eq.0) then
3457 iread=iread+lenlan+1
3458 ! print *,"iread",iread
3459 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3460 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3462 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3464 ! print *,"iread",iread
3465 if (iread.gt.lenrec) then
3470 ! print *,"ireade",ireade
3471 do while (ireade.lt.lenrec .and. &
3472 .not.iblnk(rekord(ireade:ireade)))
3475 wartosc=rekord(iread:ireade)
3477 end subroutine readsss
3478 !----------------------------------------------------------------------------
3479 subroutine multreads(rekord,lancuch,tablica,dim,default)
3482 character*(*) rekord,lancuch,tablica(dim),default
3483 character(len=80) :: aux
3484 integer :: lenlan,lenrec,iread,ireade
3491 lenlan=ilen(lancuch)
3493 iread=index(rekord,lancuch(:lenlan)//"=")
3494 ! print *,"rekord",rekord," lancuch",lancuch
3495 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3496 if (iread.eq.0) return
3497 iread=iread+lenlan+1
3499 ! print *,"iread",iread
3500 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3501 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3503 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3505 ! print *,"iread",iread
3506 if (iread.gt.lenrec) return
3508 ! print *,"ireade",ireade
3509 do while (ireade.lt.lenrec .and. &
3510 .not.iblnk(rekord(ireade:ireade)))
3513 tablica(i)=rekord(iread:ireade)
3516 end subroutine multreads
3517 !----------------------------------------------------------------------------
3518 subroutine split_string(rekord,tablica,dim,nsub)
3520 integer :: dim,nsub,i,ii,ll,kk
3521 character*(*) tablica(dim)
3522 character*(*) rekord
3532 ! Find the start of term name
3534 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
3537 ! Parse the name into TABLICA(i) until blank found
3538 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
3540 tablica(i)(kk:kk)=rekord(ii:ii)
3543 if (kk.gt.0) nsub=nsub+1
3544 if (ii.gt.ll) return
3547 end subroutine split_string
3548 !--------------------------------------------------------------------------------
3550 !--------------------------------------------------------------------------------
3551 subroutine read_compar
3553 ! Read molecular data
3555 use conform_compar, only:alloc_compar_arrays
3556 use control_data, only:pdbref
3557 use geometry_data, only:deg2rad,rad2deg
3559 ! include 'DIMENSIONS'
3560 ! include 'DIMENSIONS.ZSCOPT'
3561 ! include 'DIMENSIONS.COMPAR'
3562 ! include 'DIMENSIONS.FREE'
3563 ! include 'COMMON.IOUNITS'
3564 ! include 'COMMON.TIME1'
3565 ! include 'COMMON.SBRIDGE'
3566 ! include 'COMMON.CONTROL'
3567 ! include 'COMMON.COMPAR'
3568 ! include 'COMMON.CHAIN'
3569 ! include 'COMMON.HEADER'
3570 ! include 'COMMON.GEO'
3571 ! include 'COMMON.FREE'
3572 character(len=320) :: controlcard !,ucase
3573 character(len=64) :: wfile
3577 !elwrite(iout,*)"jestesmy w read_compar"
3578 call card_concat(controlcard,.true.)
3579 pdbref=(index(controlcard,'PDBREF').gt.0)
3580 call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
3581 call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
3582 call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
3583 call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
3584 verbose = index(controlcard,"VERBOSE").gt.0
3585 lgrp=index(controlcard,"STATIN").gt.0
3586 lgrp_out=index(controlcard,"STATOUT").gt.0
3587 merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
3588 binary = index(controlcard,"BINARY").gt.0
3589 rmscut_base_up=rmscut_base_up/50
3590 rmscut_base_low=rmscut_base_low/50
3591 call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
3592 call readi(controlcard,'NLEVEL',nlevel,1)
3593 if (nlevel.lt.0) then
3595 call alloc_compar_arrays(maxfrag,1)
3598 allocate(nfrag(nlevel))
3600 ! Read the data pertaining to elementary fragments (level 1)
3601 call readi(controlcard,'NFRAG',nfrag(1),0)
3602 write(iout,*)"nfrag(1)",nfrag(1)
3603 call alloc_compar_arrays(nfrag(1),nlevel)
3605 call card_concat(controlcard,.true.)
3606 write (iout,*) controlcard(:ilen(controlcard))
3607 call readi(controlcard,'NPIECE',npiece(j,1),0)
3608 call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
3609 call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
3610 call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
3611 call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
3612 call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
3613 call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
3614 call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
3615 call readi(controlcard,'RMS',irms(j,1),0)
3616 call readi(controlcard,'LOCAL',iloc(j),1)
3617 call readi(controlcard,'ELCONT',ielecont(j,1),1)
3618 if (ielecont(j,1).eq.0) then
3619 call readi(controlcard,'SCCONT',isccont(j,1),1)
3621 ang_cut(j)=ang_cut(j)*deg2rad
3622 ang_cut1(j)=ang_cut1(j)*deg2rad
3624 call card_concat(controlcard,.true.)
3625 call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
3626 call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
3628 write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
3629 (ifrag(1,k,j),ifrag(2,k,j),&
3630 k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
3631 " ang_cut1",ang_cut1(j)*rad2deg
3632 write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
3633 write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
3634 write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
3635 " ilocal",iloc(j)," isccont",isccont(j,1)
3637 ! Read data pertaning to higher levels
3639 call card_concat(controlcard,.true.)
3640 call readi(controlcard,'NFRAG',NFRAG(i),0)
3641 write (iout,*) "i",i," nfrag",nfrag(i)
3643 call card_concat(controlcard,.true.)
3645 call readi(controlcard,'ELCONT',ielecont(j,i),0)
3646 if (ielecont(j,i).eq.0) then
3647 call readi(controlcard,'SCCONT',isccont(j,i),1)
3649 call readi(controlcard,'RMS',irms(j,i),0)
3655 call readi(controlcard,'NPIECE',npiece(j,i),0)
3656 call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
3657 call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
3658 call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
3660 call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
3661 call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
3662 write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
3663 n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
3664 " isccont",isccont(j,i)," irms",irms(j,i)
3665 write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
3666 write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
3667 write(iout,*)"nc_frac",nc_fragm(j,i),&
3668 " nc_req",nc_req_setf(j,i)
3671 if (binary) write (iout,*) "Classes written in binary format."
3674 call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
3675 call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
3676 call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
3677 call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
3678 call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
3679 call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
3680 call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
3681 call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
3682 call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
3683 call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
3684 call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
3685 call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
3686 call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
3687 call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
3688 call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
3689 call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
3690 call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
3691 call readi(controlcard,'RMS_SINGLE',irms_single,0)
3692 call readi(controlcard,'CONT_SINGLE',icont_single,1)
3693 call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
3694 call readi(controlcard,'RMS_PAIR',irms_pair,0)
3695 call readi(controlcard,'CONT_PAIR',icont_pair,1)
3696 call readi(controlcard,'SPLIT_BET',isplit_bet,0)
3697 angcut_hel=angcut_hel*deg2rad
3698 angcut1_hel=angcut1_hel*deg2rad
3699 angcut_bet=angcut_bet*deg2rad
3700 angcut1_bet=angcut1_bet*deg2rad
3701 angcut_strand=angcut_strand*deg2rad
3702 angcut1_strand=angcut1_strand*deg2rad
3703 write (iout,*) "Automatic detection of structural elements"
3704 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
3705 ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
3706 ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
3707 ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
3708 ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
3709 ' SPLIT_BET',isplit_bet
3710 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
3711 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
3712 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
3713 ' MAXANG_HEL',angcut1_hel*rad2deg
3714 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
3715 ' MAXANG_BET',angcut1_bet*rad2deg
3716 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
3717 ' MAXANG_STRAND',angcut1_strand*rad2deg
3718 write (iout,*) 'FRAC_MIN',frac_min_set
3720 end subroutine read_compar
3721 !--------------------------------------------------------------------------------
3723 !--------------------------------------------------------------------------------
3724 subroutine read_ref_structure(*)
3726 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
3729 use control_data, only:pdbref
3730 use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
3731 nstart_sup,nstart_seq,nperm,nres0
3732 use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
3733 use compare, only:seq_comp !,contact,elecont
3734 use geometry, only:chainbuild,dist
3735 use io_config, only:readpdb
3737 use conform_compar, only:contact,elecont
3739 ! include 'DIMENSIONS'
3740 ! include 'DIMENSIONS.ZSCOPT'
3741 ! include 'DIMENSIONS.COMPAR'
3742 ! include 'COMMON.IOUNITS'
3743 ! include 'COMMON.GEO'
3744 ! include 'COMMON.VAR'
3745 ! include 'COMMON.INTERACT'
3746 ! include 'COMMON.LOCAL'
3747 ! include 'COMMON.NAMES'
3748 ! include 'COMMON.CHAIN'
3749 ! include 'COMMON.FFIELD'
3750 ! include 'COMMON.SBRIDGE'
3751 ! include 'COMMON.HEADER'
3752 ! include 'COMMON.CONTROL'
3753 ! include 'COMMON.CONTACTS1'
3754 ! include 'COMMON.PEPTCONT'
3755 ! include 'COMMON.TIME1'
3756 ! include 'COMMON.COMPAR'
3757 character(len=4) :: sequence(nres)
3759 !el real(kind=8) :: x(maxvar)
3760 integer :: itype_pdb(nres,5)
3761 !el logical seq_comp
3762 integer :: i,j,k,nres_pdb,iaux,mnum
3763 real(kind=8) :: ddsc !el,dist
3764 integer :: kkk !,ilen
3768 write (iout,*) "pdbref",pdbref
3770 read(inp,'(a)') pdbfile
3771 write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
3772 pdbfile(:ilen(pdbfile))
3773 open(ipdbin,file=pdbfile,status='old',err=33)
3775 33 write (iout,'(a)') 'Error opening PDB file.'
3780 itype_pdb(i,mnum)=itype(i,mnum)
3786 iaux=itype_pdb(i,mnum)
3787 itype_pdb(i,mnum)=itype(i,mnum)
3795 if (nsup.le.(nct-nnt+1)) then
3796 do i=0,nct-nnt+1-nsup
3797 if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
3799 do j=nnt+nsup-1,nnt,-1
3801 cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
3804 do j=nnt+nsup-1,nnt,-1
3806 cref(k,j+i,kkk)=cref(k,j,kkk)
3808 write(iout,*) "J",j,"J+I",j+i
3809 phi_ref(j+i)=phi_ref(j)
3810 theta_ref(j+i)=theta_ref(j)
3811 alph_ref(j+i)=alph_ref(j)
3812 omeg_ref(j+i)=omeg_ref(j)
3816 write (iout,'(i5,3f10.5,5x,3f10.5)') &
3817 j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
3825 write (iout,'(a)') &
3826 'Error - sequences to be superposed do not match.'
3829 do i=0,nsup-(nct-nnt+1)
3830 if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
3833 nstart_sup=nstart_sup+i
3838 write (iout,'(a)') &
3839 'Error - sequences to be superposed do not match.'
3843 write (iout,'(a,i5)') &
3844 'Experimental structure begins at residue',nstart_seq
3846 call read_angles(inp,*38)
3848 38 write (iout,'(a)') 'Error reading reference structure.'
3857 cref(j,i,kkk)=c(j,i)
3861 nend_sup=nstart_sup+nsup-1
3864 c(j,i)=cref(j,i,kkk)
3870 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
3872 if (itype(i,mnum).ne.10) then
3873 ddsc = dist(i,nres+i)
3875 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
3879 dc_norm(j,nres+i)=0.0d0
3882 ! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
3883 ! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
3884 ! dc_norm(3,nres+i)**2
3886 dc(j,i)=c(j,i+1)-c(j,i)
3890 dc_norm(j,i)=dc(j,i)/ddsc
3893 ! print *,"Calling contact"
3894 call contact(.true.,ncont_ref,icont_ref(1,1),&
3895 nstart_sup,nend_sup)
3896 ! print *,"Calling elecont"
3897 call elecont(.true.,ncont_pept_ref,&
3898 icont_pept_ref(1,1),&
3899 nstart_sup,nend_sup)
3900 write (iout,'(a,i3,a,i3,a,i3,a)') &
3901 'Number of residues to be superposed:',nsup,&
3902 ' (from residue',nstart_sup,' to residue',&
3905 end subroutine read_ref_structure
3906 !--------------------------------------------------------------------------------
3908 !--------------------------------------------------------------------------------
3909 subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
3911 use geometry_data, only:nres,c,boxxsize,boxysize,boxzsize
3912 use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
3913 use energy, only:boxshift
3914 ! implicit real*8 (a-h,o-z)
3915 ! include 'DIMENSIONS'
3916 ! include 'DIMENSIONS.ZSCOPT'
3917 ! include 'COMMON.CHAIN'
3918 ! include 'COMMON.INTERACT'
3919 ! include 'COMMON.NAMES'
3920 ! include 'COMMON.IOUNITS'
3921 ! include 'COMMON.HEADER'
3922 ! include 'COMMON.SBRIDGE'
3923 character(len=50) :: tytul
3924 character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',&
3925 'D','E','F','G','H','I','J','K','L','M','N','O',&
3926 'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
3927 integer,dimension(nres) :: ica !(maxres)
3928 real(kind=8) :: temp,efree,etot,entropy,rmsdev,xj,yj,zj
3929 integer :: ii,i,j,iti,ires,iatom,ichain,mnum
3930 write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
3932 write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
3934 write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
3942 if (iti.eq.ntyp1) then
3943 if (itype(i-1,molnum(i-1)).eq.ntyp1) then
3946 write (ipdb,'(a)') 'TER'
3953 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
3956 xj=boxshift(c(1,i)-c(1,2),boxxsize)
3957 yj=boxshift(c(2,i)-c(2,2),boxysize)
3958 zj=boxshift(c(3,i)-c(3,2),boxzsize)
3959 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
3960 ires,c(1,2)+xj,c(2,2)+yj,c(3,2)+zj
3962 if ((iti.ne.10).and.(mnum.ne.5)) then
3964 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
3965 ires,(c(j,nres+i),j=1,3)
3969 write (ipdb,'(a)') 'TER'
3972 if (itype(i,mnum).eq.ntyp1) cycle
3973 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3974 write (ipdb,30) ica(i),ica(i+1)
3975 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3976 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
3977 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
3978 write (ipdb,30) ica(i),ica(i)+1
3981 if (itype(nct,molnum(nct)).ne.10) then
3982 write (ipdb,30) ica(nct),ica(nct)+1
3985 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
3987 write (ipdb,'(a6)') 'ENDMDL'
3988 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3989 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3990 30 FORMAT ('CONECT',8I5)
3992 end subroutine pdboutW
3994 !------------------------------------------------------------------------------
3996 !-----------------------------------------------------------------------------
3997 !-----------------------------------------------------------------------------