Adam's changes
[unres.git] / source / wham / src-restraints / parmread.F
1       subroutine parmread(iparm,*)
2 C
3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
5 C
6       implicit real*8 (a-h,o-z)
7       include 'DIMENSIONS'
8       include 'DIMENSIONS.ZSCOPT'
9       include 'DIMENSIONS.FREE'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.CHAIN'
12       include 'COMMON.INTERACT'
13       include 'COMMON.GEO'
14       include 'COMMON.LOCAL'
15       include 'COMMON.TORSION'
16       include 'COMMON.FFIELD'
17       include 'COMMON.NAMES'
18       include 'COMMON.SBRIDGE'
19       include 'COMMON.WEIGHTS'
20       include 'COMMON.ENEPS'
21       include 'COMMON.SCCOR'
22       include 'COMMON.SCROT'
23       include 'COMMON.FREE'
24       character*1 t1,t2,t3
25       character*1 onelett(4) /"G","A","P","D"/
26       logical lprint
27       dimension blower(3,3,maxlob)
28       character*800 controlcard
29       character*256 bondname_t,thetname_t,rotname_t,torname_t,
30      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
31      & sccorname_t
32       integer ilen
33       external ilen
34       character*16 key
35       integer iparm
36       double precision ip,mp
37 C
38 C Body
39 C
40 C Set LPRINT=.TRUE. for debugging
41       dwa16=2.0d0**(1.0d0/6.0d0)
42       lprint=.false.
43       itypro=20
44 C Assign virtual-bond length
45       vbl=3.8D0
46       vblinv=1.0D0/vbl
47       vblinv2=vblinv*vblinv
48       call card_concat(controlcard,.true.)
49       wname(4)="WCORRH"
50       do i=1,n_ene
51         key = wname(i)(:ilen(wname(i)))
52         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
53       enddo
54       call reada(controlcard,"D0CM",d0cm,3.78d0)
55       call reada(controlcard,"AKCM",akcm,15.1d0)
56       call reada(controlcard,"AKTH",akth,11.0d0)
57       call reada(controlcard,"AKCT",akct,12.0d0)
58       call reada(controlcard,"V1SS",v1ss,-1.08d0)
59       call reada(controlcard,"V2SS",v2ss,7.61d0)
60       call reada(controlcard,"V3SS",v3ss,13.7d0)
61       call reada(controlcard,"EBR",ebr,-5.50D0)
62 c      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
63       write (iout,*) "iparm",iparm," myparm",myparm
64 c      do i=1,maxres
65 c        dyn_ss_mask(i)=.false.
66 c      enddo
67       do i=1,maxres-1
68         do j=i+1,maxres
69           dyn_ssbond_ij(i,j)=1.0d300
70         enddo
71       enddo
72       call reada(controlcard,"HT",Ht,0.0D0)
73
74 c      if(me.eq.king.or..not.out1file) then
75 c       print *,'indpdb=',indpdb,' pdbref=',pdbref
76 c      endif
77 c If reading not own parameters, skip assignment
78 cc      write(iout,*) "KURWA", ww(15)
79       if (iparm.eq.myparm .or. .not.separate_parset) then
80
81 c
82 c Setup weights for UNRES
83 c
84       wsc=ww(1)
85       wscp=ww(2)
86       welec=ww(3)
87       wcorr=ww(4)
88       wcorr5=ww(5)
89       wcorr6=ww(6)
90       wel_loc=ww(7)
91       wturn3=ww(8)
92       wturn4=ww(9)
93       wturn6=ww(10)
94       wang=ww(11)
95       wscloc=ww(12)
96       wtor=ww(13)
97       wtor_d=ww(14)
98       wvdwpp=ww(16)
99       wstrain=ww(15)
100       wbond=ww(18)
101       wsccor=ww(19)
102
103       endif
104 cc      write(iout,*) "KURWA", wstrain,akcm,akth,wsc,dyn_ss
105
106       call card_concat(controlcard,.false.)
107
108 c Return if not own parameters
109
110       if (iparm.ne.myparm .and. separate_parset) return
111
112       call reads(controlcard,"BONDPAR",bondname_t,bondname)
113       open (ibond,file=bondname_t,status='old')
114       rewind(ibond)
115       call reads(controlcard,"THETPAR",thetname_t,thetname)
116       open (ithep,file=thetname_t,status='old')
117       rewind(ithep) 
118       call reads(controlcard,"ROTPAR",rotname_t,rotname)
119       open (irotam,file=rotname_t,status='old')
120       rewind(irotam)
121       call reads(controlcard,"TORPAR",torname_t,torname)
122       open (itorp,file=torname_t,status='old')
123       rewind(itorp)
124       call reads(controlcard,"TORDPAR",tordname_t,tordname)
125       open (itordp,file=tordname_t,status='old')
126       rewind(itordp)
127       call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
128       open (isccor,file=sccorname_t,status='old')
129       rewind(isccor)
130       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
131       open (ifourier,file=fouriername_t,status='old')
132       rewind(ifourier)
133       call reads(controlcard,"ELEPAR",elename_t,elename)
134       open (ielep,file=elename_t,status='old')
135       rewind(ielep)
136       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
137       open (isidep,file=sidename_t,status='old')
138       rewind(isidep)
139       call reads(controlcard,"SCPPAR",scpname_t,scpname)
140       open (iscpp,file=scpname_t,status='old')
141       rewind(iscpp)
142       write (iout,*) "Parameter set:",iparm
143       write (iout,*) "Energy-term weights:"
144       do i=1,n_ene
145         write (iout,'(a16,f10.5)') wname(i),ww(i)
146       enddo
147       write (iout,*) "Sidechain potential file        : ",
148      &  sidename_t(:ilen(sidename_t))
149 #ifndef OLDSCP
150       write (iout,*) "SCp potential file              : ",
151      &  scpname_t(:ilen(scpname_t))
152 #endif  
153       write (iout,*) "Electrostatic potential file    : ",
154      &  elename_t(:ilen(elename_t))
155       write (iout,*) "Cumulant coefficient file       : ",
156      &  fouriername_t(:ilen(fouriername_t))
157       write (iout,*) "Torsional parameter file        : ",
158      &  torname_t(:ilen(torname_t))
159       write (iout,*) "Double torsional parameter file : ",
160      &  tordname_t(:ilen(tordname_t))
161       write (iout,*) "Backbone-rotamer parameter file : ",
162      &  sccorname(:ilen(sccorname))
163       write (iout,*) "Bond & inertia constant file    : ",
164      &  bondname_t(:ilen(bondname_t))
165       write (iout,*) "Bending parameter file          : ",
166      &  thetname_t(:ilen(thetname_t))
167       write (iout,*) "Rotamer parameter file          : ",
168      &  rotname_t(:ilen(rotname_t))
169
170 c
171 c Read the virtual-bond parameters, masses, and moments of inertia
172 c and Stokes' radii of the peptide group and side chains
173 c
174 #ifdef CRYST_BOND
175       read (ibond,*) vbldp0,akp
176       do i=1,ntyp
177         nbondterm(i)=1
178         read (ibond,*) vbldsc0(1,i),aksc(1,i)
179         dsc(i) = vbldsc0(1,i)
180         if (i.eq.10) then
181           dsc_inv(i)=0.0D0
182         else
183           dsc_inv(i)=1.0D0/dsc(i)
184         endif
185       enddo
186 #else
187       read (ibond,*) ijunk,vbldp0,akp,rjunk
188       do i=1,ntyp
189         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
190      &   j=1,nbondterm(i))
191         dsc(i) = vbldsc0(1,i)
192         if (i.eq.10) then
193           dsc_inv(i)=0.0D0
194         else
195           dsc_inv(i)=1.0D0/dsc(i)
196         endif
197       enddo
198 #endif
199       if (lprint) then
200         write(iout,'(/a/)')"Force constants virtual bonds:"
201         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
202      &   'inertia','Pstok'
203         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
204         do i=1,ntyp
205           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
206      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
207           do j=2,nbondterm(i)
208             write (iout,'(13x,3f10.5)')
209      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
210           enddo
211         enddo
212       endif
213 #ifdef CRYST_THETA
214 C
215 C Read the parameters of the probability distribution/energy expression 
216 C of the virtual-bond valence angles theta
217 C
218       do i=1,ntyp
219         read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
220         read (ithep,*) (polthet(j,i),j=0,3)
221         read (ithep,*) (gthet(j,i),j=1,3)
222         read (ithep,*) theta0(i),sig0(i),sigc0(i)
223         sigc0(i)=sigc0(i)**2
224       enddo
225       close (ithep)
226       if (lprint) then
227 c       write (iout,'(a)') 
228 c    &    'Parameters of the virtual-bond valence angles:'
229 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
230 c    & '    ATHETA0   ','         A1   ','        A2    ',
231 c    & '        B1    ','         B2   '        
232 c       do i=1,ntyp
233 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
234 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
235 c       enddo
236 c       write (iout,'(/a/9x,5a/79(1h-))') 
237 c    & 'Parameters of the expression for sigma(theta_c):',
238 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
239 c    & '     ALPH3    ','    SIGMA0C   '        
240 c       do i=1,ntyp
241 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
242 c    &      (polthet(j,i),j=0,3),sigc0(i) 
243 c       enddo
244 c       write (iout,'(/a/9x,5a/79(1h-))') 
245 c    & 'Parameters of the second gaussian:',
246 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
247 c    & '        G2    ','         G3   '        
248 c       do i=1,ntyp
249 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
250 c    &       sig0(i),(gthet(j,i),j=1,3)
251 c       enddo
252         write (iout,'(a)') 
253      &    'Parameters of the virtual-bond valence angles:'
254         write (iout,'(/a/9x,5a/79(1h-))') 
255      & 'Coefficients of expansion',
256      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
257      & '   b1*10^1    ','    b2*10^1   '        
258         do i=1,ntyp
259           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
260      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
261         enddo
262         write (iout,'(/a/9x,5a/79(1h-))') 
263      & 'Parameters of the expression for sigma(theta_c):',
264      & ' alpha0       ','  alph1       ',' alph2        ',
265      & ' alhp3        ','   sigma0c    '        
266         do i=1,ntyp
267           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
268      &      (polthet(j,i),j=0,3),sigc0(i) 
269         enddo
270         write (iout,'(/a/9x,5a/79(1h-))') 
271      & 'Parameters of the second gaussian:',
272      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
273      & '        G2    ','   G3*10^1    '        
274         do i=1,ntyp
275           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
276      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
277         enddo
278       endif
279 #else
280 C
281 C Read the parameters of Utheta determined from ab initio surfaces
282 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
283 C
284       read (ithep,*) nthetyp,ntheterm,ntheterm2,
285      &  ntheterm3,nsingle,ndouble
286       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
287       read (ithep,*) (ithetyp(i),i=1,ntyp1)
288       do i=1,maxthetyp
289         do j=1,maxthetyp
290           do k=1,maxthetyp
291             aa0thet(i,j,k)=0.0d0
292             do l=1,ntheterm
293               aathet(l,i,j,k)=0.0d0
294             enddo
295             do l=1,ntheterm2
296               do m=1,nsingle
297                 bbthet(m,l,i,j,k)=0.0d0
298                 ccthet(m,l,i,j,k)=0.0d0
299                 ddthet(m,l,i,j,k)=0.0d0
300                 eethet(m,l,i,j,k)=0.0d0
301               enddo
302             enddo
303             do l=1,ntheterm3
304               do m=1,ndouble
305                 do mm=1,ndouble
306                  ffthet(mm,m,l,i,j,k)=0.0d0
307                  ggthet(mm,m,l,i,j,k)=0.0d0
308                 enddo
309               enddo
310             enddo
311           enddo
312         enddo
313       enddo
314       do i=1,nthetyp
315         do j=1,nthetyp
316           do k=1,nthetyp
317             read (ithep,'(3a)') res1,res2,res3
318             read (ithep,*) aa0thet(i,j,k)
319             read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
320             read (ithep,*)
321      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
322      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
323      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
324      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
325             read (ithep,*)
326      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
327      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
328      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
329           enddo
330         enddo
331       enddo
332 C
333 C For dummy ends assign glycine-type coefficients of theta-only terms; the
334 C coefficients of theta-and-gamma-dependent terms are zero.
335 C
336       do i=1,nthetyp
337         do j=1,nthetyp
338           do l=1,ntheterm
339             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
340             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
341           enddo
342           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
343           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
344         enddo
345         do l=1,ntheterm
346           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
347         enddo
348         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
349       enddo
350 C
351 C Control printout of the coefficients of virtual-bond-angle potentials
352 C
353       if (lprint) then
354         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
355         do i=1,nthetyp+1
356           do j=1,nthetyp+1
357             do k=1,nthetyp+1
358               write (iout,'(//4a)')
359      &         'Type ',onelett(i),onelett(j),onelett(k)
360               write (iout,'(//a,10x,a)') " l","a[l]"
361               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
362               write (iout,'(i2,1pe15.5)')
363      &           (l,aathet(l,i,j,k),l=1,ntheterm)
364             do l=1,ntheterm2
365               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
366      &          "b",l,"c",l,"d",l,"e",l
367               do m=1,nsingle
368                 write (iout,'(i2,4(1pe15.5))') m,
369      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
370      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
371               enddo
372             enddo
373             do l=1,ntheterm3
374               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
375      &          "f+",l,"f-",l,"g+",l,"g-",l
376               do m=2,ndouble
377                 do n=1,m-1
378                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
379      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
380      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
381                 enddo
382               enddo
383             enddo
384           enddo
385         enddo
386       enddo
387       call flush(iout)
388       endif
389 #endif
390
391 #ifdef CRYST_SC
392 C
393 C Read the parameters of the probability distribution/energy expression
394 C of the side chains.
395 C
396       do i=1,ntyp
397         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
398         if (i.eq.10) then
399           dsc_inv(i)=0.0D0
400         else
401           dsc_inv(i)=1.0D0/dsc(i)
402         endif
403         if (i.ne.10) then
404         do j=1,nlob(i)
405           do k=1,3
406             do l=1,3
407               blower(l,k,j)=0.0D0
408             enddo
409           enddo
410         enddo  
411         bsc(1,i)=0.0D0
412         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
413         do j=2,nlob(i)
414           read (irotam,*) bsc(j,i)
415           read (irotam,*) (censc(k,j,i),k=1,3),
416      &                                 ((blower(k,l,j),l=1,k),k=1,3)
417         enddo
418         do j=1,nlob(i)
419           do k=1,3
420             do l=1,k
421               akl=0.0D0
422               do m=1,3
423                 akl=akl+blower(k,m,j)*blower(l,m,j)
424               enddo
425               gaussc(k,l,j,i)=akl
426               gaussc(l,k,j,i)=akl
427             enddo
428           enddo 
429         enddo
430         endif
431       enddo
432       close (irotam)
433       if (lprint) then
434         write (iout,'(/a)') 'Parameters of side-chain local geometry'
435         do i=1,ntyp
436           nlobi=nlob(i)
437           if (nlobi.gt.0) then
438           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
439      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
440 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
441 c          write (iout,'(a,f10.4,4(16x,f10.4))')
442 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
443 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
444            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
445      &                             'log h',(bsc(j,i),j=1,nlobi)
446            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
447      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
448 c          write (iout,'(a)')
449 c         do j=1,nlobi
450 c           ind=0
451 c           do k=1,3
452 c             do l=1,k
453 c              ind=ind+1
454 c              blower(k,l,j)=gaussc(ind,j,i)
455 c             enddo
456 c           enddo
457 c         enddo
458           do k=1,3
459             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
460      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
461           enddo
462           endif
463         enddo
464       endif
465 #else
466 C
467 C Read scrot parameters for potentials determined from all-atom AM1 calculations
468 C added by Urszula Kozlowska 07/11/2007
469 C
470       do i=1,ntyp
471         read (irotam,*)
472        if (i.eq.10) then
473          read (irotam,*)
474        else
475          do j=1,65
476            read(irotam,*) sc_parmin(j,i)
477          enddo
478        endif
479       enddo
480 #endif
481       close(irotam)
482 #ifdef CRYST_TOR
483 C
484 C Read torsional parameters in old format
485 C
486       read (itorp,*) ntortyp,nterm_old
487       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
488       read (itorp,*) (itortyp(i),i=1,ntyp)
489       do i=1,ntortyp
490         do j=1,ntortyp
491           read (itorp,'(a)')
492           do k=1,nterm_old
493             read (itorp,*) kk,v1(k,j,i),v2(k,j,i) 
494           enddo
495         enddo
496       enddo
497       close (itorp)
498       if (lprint) then
499         write (iout,'(/a/)') 'Torsional constants:'
500         do i=1,ntortyp
501           do j=1,ntortyp
502             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
503             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
504           enddo
505         enddo
506       endif
507
508
509 #else
510 C
511 C Read torsional parameters
512 C
513       read (itorp,*) ntortyp
514       read (itorp,*) (itortyp(i),i=1,ntyp)
515       write (iout,*) 'ntortyp',ntortyp
516       do i=1,ntortyp
517         do j=1,ntortyp
518           read (itorp,*) nterm(i,j),nlor(i,j)
519           v0ij=0.0d0
520           si=-1.0d0
521           do k=1,nterm(i,j)
522             read (itorp,*) kk,v1(k,i,j),v2(k,i,j) 
523             v0ij=v0ij+si*v1(k,i,j)
524             si=-si
525           enddo
526           do k=1,nlor(i,j)
527             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) 
528             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
529           enddo
530           v0(i,j)=v0ij
531         enddo
532       enddo
533       close (itorp)
534       if (lprint) then
535         write (iout,'(/a/)') 'Torsional constants:'
536         do i=1,ntortyp
537           do j=1,ntortyp
538             write (iout,*) 'ityp',i,' jtyp',j
539             write (iout,*) 'Fourier constants'
540             do k=1,nterm(i,j)
541               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
542             enddo
543             write (iout,*) 'Lorenz constants'
544             do k=1,nlor(i,j)
545               write (iout,'(3(1pe15.5))') 
546      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
547             enddo
548           enddo
549         enddo
550       endif
551 C
552 C 6/23/01 Read parameters for double torsionals
553 C
554       do i=1,ntortyp
555         do j=1,ntortyp
556           do k=1,ntortyp
557             read (itordp,'(3a1)') t1,t2,t3
558             if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
559      &        .or. t3.ne.onelett(k)) then
560               write (iout,*) "Error in double torsional parameter file",
561      &         i,j,k,t1,t2,t3
562                stop "Error in double torsional parameter file"
563             endif
564             read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
565             read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
566             read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
567             read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
568             read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
569             read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
570      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
571           enddo
572         enddo
573       enddo
574       if (lprint) then
575       write (iout,*) 
576       write (iout,*) 'Constants for double torsionals'
577       do i=1,ntortyp
578         do j=1,ntortyp 
579           do k=1,ntortyp
580             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
581      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
582             write (iout,*)
583             write (iout,*) 'Single angles:'
584             do l=1,ntermd_1(i,j,k)
585               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
586      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
587      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
588             enddo
589             write (iout,*)
590             write (iout,*) 'Pairs of angles:'
591             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
592             do l=1,ntermd_2(i,j,k)
593               write (iout,'(i5,20f10.5)') 
594      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
595             enddo
596             write (iout,*)
597             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
598             do l=1,ntermd_2(i,j,k)
599               write (iout,'(i5,20f10.5)') 
600      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
601             enddo
602             write (iout,*)
603           enddo
604         enddo
605       enddo
606       endif
607 #endif
608 C Read of Side-chain backbone correlation parameters
609 C Modified 11 May 2012 by Adasko
610 CCC
611 C
612       read (isccor,*) nsccortyp
613       read (isccor,*) (isccortyp(i),i=1,ntyp)
614 c      write (iout,*) 'ntortyp',ntortyp
615       maxinter=3
616 cc maxinter is maximum interaction sites
617       do l=1,maxinter
618       do i=1,nsccortyp
619         do j=1,nsccortyp
620           read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
621           v0ijsccor=0.0d0
622           si=-1.0d0
623
624           do k=1,nterm_sccor(i,j)
625             read (isccor,*) kk,v1sccor(k,l,i,j)
626      &    ,v2sccor(k,l,i,j)
627             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
628             si=-si
629           enddo
630           do k=1,nlor_sccor(i,j)
631             read (isccor,*) kk,vlor1sccor(k,i,j),
632      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
633             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
634      &(1+vlor3sccor(k,i,j)**2)
635           enddo
636           v0sccor(i,j)=v0ijsccor
637         enddo
638       enddo
639       enddo
640       close (isccor)
641
642       if (lprint) then
643         write (iout,'(/a/)') 'Torsional constants:'
644         do i=1,nsccortyp
645           do j=1,nsccortyp
646             write (iout,*) 'ityp',i,' jtyp',j
647             write (iout,*) 'Fourier constants'
648             do k=1,nterm_sccor(i,j)
649            write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
650             enddo
651             write (iout,*) 'Lorenz constants'
652             do k=1,nlor_sccor(i,j)
653               write (iout,'(3(1pe15.5))')
654      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
655             enddo
656           enddo
657         enddo
658       endif
659
660 C
661 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
662 C         interaction energy of the Gly, Ala, and Pro prototypes.
663 C
664       read (ifourier,*) nloctyp
665       do i=1,nloctyp
666         read (ifourier,*)
667         read (ifourier,*) (b(ii,i),ii=1,13)
668         if (lprint) then
669         write (iout,*) 'Type',i
670         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
671         endif
672         B1(1,i)  = b(3,i)
673         B1(2,i)  = b(5,i)
674         B1tilde(1,i) = b(3,i)
675         B1tilde(2,i) =-b(5,i) 
676         B2(1,i)  = b(2,i)
677         B2(2,i)  = b(4,i)
678         CC(1,1,i)= b(7,i)
679         CC(2,2,i)=-b(7,i)
680         CC(2,1,i)= b(9,i)
681         CC(1,2,i)= b(9,i)
682         Ctilde(1,1,i)=b(7,i)
683         Ctilde(1,2,i)=b(9,i)
684         Ctilde(2,1,i)=-b(9,i)
685         Ctilde(2,2,i)=b(7,i)
686         DD(1,1,i)= b(6,i)
687         DD(2,2,i)=-b(6,i)
688         DD(2,1,i)= b(8,i)
689         DD(1,2,i)= b(8,i)
690         Dtilde(1,1,i)=b(6,i)
691         Dtilde(1,2,i)=b(8,i)
692         Dtilde(2,1,i)=-b(8,i)
693         Dtilde(2,2,i)=b(6,i)
694         EE(1,1,i)= b(10,i)+b(11,i)
695         EE(2,2,i)=-b(10,i)+b(11,i)
696         EE(2,1,i)= b(12,i)-b(13,i)
697         EE(1,2,i)= b(12,i)+b(13,i)
698       enddo
699       if (lprint) then
700       do i=1,nloctyp
701         write (iout,*) 'Type',i
702         write (iout,*) 'B1'
703 c        write (iout,'(f10.5)') B1(:,i)
704         write(iout,*) B1(1,i),B1(2,i)
705         write (iout,*) 'B2'
706 c        write (iout,'(f10.5)') B2(:,i)
707         write(iout,*) B2(1,i),B2(2,i)
708         write (iout,*) 'CC'
709         do j=1,2
710           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
711         enddo
712         write(iout,*) 'DD'
713         do j=1,2
714           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
715         enddo
716         write(iout,*) 'EE'
717         do j=1,2
718           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
719         enddo
720       enddo
721       endif
722
723 C Read electrostatic-interaction parameters
724 C
725       if (lprint) then
726         write (iout,'(/a)') 'Electrostatic interaction constants:'
727         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
728      &            'IT','JT','APP','BPP','AEL6','AEL3'
729       endif
730       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
731       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
732       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
733       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
734       close (ielep)
735       do i=1,2
736         do j=1,2
737         rri=rpp(i,j)**6
738         app (i,j)=epp(i,j)*rri*rri 
739         bpp (i,j)=-2.0D0*epp(i,j)*rri
740         ael6(i,j)=elpp6(i,j)*4.2D0**6
741         ael3(i,j)=elpp3(i,j)*4.2D0**3
742         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
743      &                    ael6(i,j),ael3(i,j)
744         enddo
745       enddo
746 C
747 C Read side-chain interaction parameters.
748 C
749       read (isidep,*) ipot,expon
750       if (ipot.lt.1 .or. ipot.gt.5) then
751         write (iout,'(2a)') 'Error while reading SC interaction',
752      &               'potential file - unknown potential type.'
753         stop
754       endif
755       expon2=expon/2
756       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
757      & ', exponents are ',expon,2*expon 
758       goto (10,20,30,30,40) ipot
759 C----------------------- LJ potential ---------------------------------
760    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
761       if (lprint) then
762         write (iout,'(/a/)') 'Parameters of the LJ potential:'
763         write (iout,'(a/)') 'The epsilon array:'
764         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
765         write (iout,'(/a)') 'One-body parameters:'
766         write (iout,'(a,4x,a)') 'residue','sigma'
767         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
768       endif
769       goto 50
770 C----------------------- LJK potential --------------------------------
771    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
772      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
773       if (lprint) then
774         write (iout,'(/a/)') 'Parameters of the LJK potential:'
775         write (iout,'(a/)') 'The epsilon array:'
776         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
777         write (iout,'(/a)') 'One-body parameters:'
778         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
779         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
780      &        i=1,ntyp)
781       endif
782       goto 50
783 C---------------------- GB or BP potential -----------------------------
784    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
785      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
786      &  (alp(i),i=1,ntyp)
787 C For the GB potential convert sigma'**2 into chi'
788       if (ipot.eq.4) then
789         do i=1,ntyp
790           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
791         enddo
792       endif
793       if (lprint) then
794         write (iout,'(/a/)') 'Parameters of the BP potential:'
795         write (iout,'(a/)') 'The epsilon array:'
796         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
797         write (iout,'(/a)') 'One-body parameters:'
798         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
799      &       '    chip  ','    alph  '
800         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
801      &                     chip(i),alp(i),i=1,ntyp)
802       endif
803       goto 50
804 C--------------------- GBV potential -----------------------------------
805    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
806      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
807      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
808       if (lprint) then
809         write (iout,'(/a/)') 'Parameters of the GBV potential:'
810         write (iout,'(a/)') 'The epsilon array:'
811         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
812         write (iout,'(/a)') 'One-body parameters:'
813         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
814      &      's||/s_|_^2','    chip  ','    alph  '
815         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
816      &           sigii(i),chip(i),alp(i),i=1,ntyp)
817       endif
818    50 continue
819       close (isidep)
820 C-----------------------------------------------------------------------
821 C Calculate the "working" parameters of SC interactions.
822       do i=2,ntyp
823         do j=1,i-1
824           eps(i,j)=eps(j,i)
825         enddo
826       enddo
827       do i=1,ntyp
828         do j=i,ntyp
829           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
830           sigma(j,i)=sigma(i,j)
831           rs0(i,j)=dwa16*sigma(i,j)
832           rs0(j,i)=rs0(i,j)
833         enddo
834       enddo
835       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
836      & 'Working parameters of the SC interactions:',
837      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
838      & '  chi1   ','   chi2   ' 
839       do i=1,ntyp
840         do j=i,ntyp
841           epsij=eps(i,j)
842           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
843             rrij=sigma(i,j)
844           else
845             rrij=rr0(i)+rr0(j)
846           endif
847           r0(i,j)=rrij
848           r0(j,i)=rrij
849           rrij=rrij**expon
850           epsij=eps(i,j)
851           sigeps=dsign(1.0D0,epsij)
852           epsij=dabs(epsij)
853           aa(i,j)=epsij*rrij*rrij
854           bb(i,j)=-sigeps*epsij*rrij
855           aa(j,i)=aa(i,j)
856           bb(j,i)=bb(i,j)
857           if (ipot.gt.2) then
858             sigt1sq=sigma0(i)**2
859             sigt2sq=sigma0(j)**2
860             sigii1=sigii(i)
861             sigii2=sigii(j)
862             ratsig1=sigt2sq/sigt1sq
863             ratsig2=1.0D0/ratsig1
864             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
865             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
866             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
867           else
868             rsum_max=sigma(i,j)
869           endif
870 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
871             sigmaii(i,j)=rsum_max
872             sigmaii(j,i)=rsum_max 
873 c         else
874 c           sigmaii(i,j)=r0(i,j)
875 c           sigmaii(j,i)=r0(i,j)
876 c         endif
877 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
878           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
879             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
880             augm(i,j)=epsij*r_augm**(2*expon)
881 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
882             augm(j,i)=augm(i,j)
883           else
884             augm(i,j)=0.0D0
885             augm(j,i)=0.0D0
886           endif
887           if (lprint) then
888             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
889      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
890      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
891           endif
892         enddo
893       enddo
894 C
895 C Define the SC-p interaction constants and SS bond potentials
896 C
897       if (dyn_ss) then
898         ss_depth=ebr/wsc-0.25*eps(1,1)
899         Ht=Ht/wsc-0.25*eps(1,1)
900         akcm=akcm*wstrain/wsc
901         akth=akth*wstrain/wsc
902         akct=akct*wstrain/wsc
903         v1ss=v1ss*wstrain/wsc
904         v2ss=v2ss*wstrain/wsc
905         v3ss=v3ss*wstrain/wsc
906       else
907         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
908       endif
909        write (iout,*) "Parameters of the SS-bond potential:"
910        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
911      & " AKCT",akct
912        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
913        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
914        write (iout,*)" HT",Ht
915
916 #ifdef OLDSCP
917       do i=1,20
918 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
919 C helix formation)
920 c       aad(i,1)=0.3D0*4.0D0**12
921 C Following line for constants currently implemented
922 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
923         aad(i,1)=1.5D0*4.0D0**12
924 c       aad(i,1)=0.17D0*5.6D0**12
925         aad(i,2)=aad(i,1)
926 C "Soft" SC-p repulsion
927         bad(i,1)=0.0D0
928 C Following line for constants currently implemented
929 c       aad(i,1)=0.3D0*4.0D0**6
930 C "Hard" SC-p repulsion
931         bad(i,1)=3.0D0*4.0D0**6
932 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
933         bad(i,2)=bad(i,1)
934 c       aad(i,1)=0.0D0
935 c       aad(i,2)=0.0D0
936 c       bad(i,1)=1228.8D0
937 c       bad(i,2)=1228.8D0
938       enddo
939 #else
940 C
941 C 8/9/01 Read the SC-p interaction constants from file
942 C
943       do i=1,ntyp
944         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
945       enddo
946       do i=1,ntyp
947         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
948         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
949         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
950         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
951       enddo
952
953       if (lprint) then
954         write (iout,*) "Parameters of SC-p interactions:"
955         do i=1,20
956           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
957      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
958         enddo
959       endif
960 #endif
961 C
962 C Define the constants of the disulfide bridge
963 C
964 c      ebr=-5.50D0
965 c
966 c Old arbitrary potential - commented out.
967 c
968 c      dbr= 4.20D0
969 c      fbr= 3.30D0
970 c
971 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
972 c energy surface of diethyl disulfide.
973 c A. Liwo and U. Kozlowska, 11/24/03
974 c
975 c      D0CM = 3.78d0
976 c      AKCM = 15.1d0
977 c      AKTH = 11.0d0
978 c      AKCT = 12.0d0
979 c      V1SS =-1.08d0
980 c      V2SS = 7.61d0
981 c      V3SS = 13.7d0
982
983 c      if (lprint) then
984 c      write (iout,'(/a)') "Disulfide bridge parameters:"
985 c      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
986 c      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
987 c      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
988 c      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
989 c     & ' v3ss:',v3ss
990 c      endif
991       return
992       end