8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
870 if (oldion.eq.1) then
872 read(iion,*) msc(i,5),restok(i,5)
873 print *,msc(i,5),restok(i,5)
877 read (iion,*) ncatprotparm
878 allocate(catprm(ncatprotparm,4))
880 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
884 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
885 !----------------------------------------------------
886 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
887 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
888 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
889 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
890 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
891 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
901 allocate(liptranene(ntyp))
902 !C reading lipid parameters
903 write (iout,*) "iliptranpar",iliptranpar
905 read(iliptranpar,*) pepliptran
908 read(iliptranpar,*) liptranene(i)
909 print *,liptranene(i)
915 ! Read the parameters of the probability distribution/energy expression
916 ! of the virtual-bond valence angles theta
919 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
920 (bthet(j,i,1,1),j=1,2)
921 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
922 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
923 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
927 athet(1,i,1,-1)=athet(1,i,1,1)
928 athet(2,i,1,-1)=athet(2,i,1,1)
929 bthet(1,i,1,-1)=-bthet(1,i,1,1)
930 bthet(2,i,1,-1)=-bthet(2,i,1,1)
931 athet(1,i,-1,1)=-athet(1,i,1,1)
932 athet(2,i,-1,1)=-athet(2,i,1,1)
933 bthet(1,i,-1,1)=bthet(1,i,1,1)
934 bthet(2,i,-1,1)=bthet(2,i,1,1)
938 athet(1,i,-1,-1)=athet(1,-i,1,1)
939 athet(2,i,-1,-1)=-athet(2,-i,1,1)
940 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
941 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
942 athet(1,i,-1,1)=athet(1,-i,1,1)
943 athet(2,i,-1,1)=-athet(2,-i,1,1)
944 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
945 bthet(2,i,-1,1)=bthet(2,-i,1,1)
946 athet(1,i,1,-1)=-athet(1,-i,1,1)
947 athet(2,i,1,-1)=athet(2,-i,1,1)
948 bthet(1,i,1,-1)=bthet(1,-i,1,1)
949 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
954 polthet(j,i)=polthet(j,-i)
957 gthet(j,i)=gthet(j,-i)
965 'Parameters of the virtual-bond valence angles:'
966 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
967 ' ATHETA0 ',' A1 ',' A2 ',&
970 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
971 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
973 write (iout,'(/a/9x,5a/79(1h-))') &
974 'Parameters of the expression for sigma(theta_c):',&
975 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
976 ' ALPH3 ',' SIGMA0C '
978 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
979 (polthet(j,i),j=0,3),sigc0(i)
981 write (iout,'(/a/9x,5a/79(1h-))') &
982 'Parameters of the second gaussian:',&
983 ' THETA0 ',' SIGMA0 ',' G1 ',&
986 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
987 sig0(i),(gthet(j,i),j=1,3)
991 'Parameters of the virtual-bond valence angles:'
992 write (iout,'(/a/9x,5a/79(1h-))') &
993 'Coefficients of expansion',&
994 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
995 ' b1*10^1 ',' b2*10^1 '
997 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
998 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
999 (10*bthet(j,i,1,1),j=1,2)
1001 write (iout,'(/a/9x,5a/79(1h-))') &
1002 'Parameters of the expression for sigma(theta_c):',&
1003 ' alpha0 ',' alph1 ',' alph2 ',&
1004 ' alhp3 ',' sigma0c '
1006 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1007 (polthet(j,i),j=0,3),sigc0(i)
1009 write (iout,'(/a/9x,5a/79(1h-))') &
1010 'Parameters of the second gaussian:',&
1011 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1014 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1015 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1021 ! Read the parameters of Utheta determined from ab initio surfaces
1022 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1024 IF (tor_mode.eq.0) THEN
1025 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1026 ntheterm3,nsingle,ndouble
1027 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1029 !----------------------------------------------------
1030 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1031 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1032 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1033 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1034 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1035 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1036 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1037 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1038 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1039 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1040 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1041 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1042 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1043 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1044 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1045 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1046 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1047 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1048 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1049 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1050 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1051 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1052 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1053 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1055 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1057 ithetyp(i)=-ithetyp(-i)
1060 aa0thet(:,:,:,:)=0.0d0
1061 aathet(:,:,:,:,:)=0.0d0
1062 bbthet(:,:,:,:,:,:)=0.0d0
1063 ccthet(:,:,:,:,:,:)=0.0d0
1064 ddthet(:,:,:,:,:,:)=0.0d0
1065 eethet(:,:,:,:,:,:)=0.0d0
1066 ffthet(:,:,:,:,:,:,:)=0.0d0
1067 ggthet(:,:,:,:,:,:,:)=0.0d0
1069 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1071 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1072 ! VAR:1=non-glicyne non-proline 2=proline
1073 ! VAR:negative values for D-aminoacid
1075 do j=-nthetyp,nthetyp
1076 do k=-nthetyp,nthetyp
1077 read (ithep,'(6a)',end=111,err=111) res1
1078 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1079 ! VAR: aa0thet is variable describing the average value of Foureir
1080 ! VAR: expansion series
1081 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1082 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1083 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1084 read (ithep,*,end=111,err=111) &
1085 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1086 read (ithep,*,end=111,err=111) &
1087 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1088 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1089 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1090 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1092 read (ithep,*,end=111,err=111) &
1093 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1094 ffthet(lll,llll,ll,i,j,k,iblock),&
1095 ggthet(llll,lll,ll,i,j,k,iblock),&
1096 ggthet(lll,llll,ll,i,j,k,iblock),&
1097 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1102 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1103 ! coefficients of theta-and-gamma-dependent terms are zero.
1104 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1105 ! RECOMENTDED AFTER VERSION 3.3)
1109 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1110 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1112 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1113 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1116 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1118 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1121 ! AND COMMENT THE LOOPS BELOW
1125 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1126 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1128 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1129 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1132 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1134 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1139 ! Substitution for D aminoacids from symmetry.
1142 do j=-nthetyp,nthetyp
1143 do k=-nthetyp,nthetyp
1144 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1146 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1150 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1151 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1152 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1153 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1159 ffthet(llll,lll,ll,i,j,k,iblock)= &
1160 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1161 ffthet(lll,llll,ll,i,j,k,iblock)= &
1162 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1163 ggthet(llll,lll,ll,i,j,k,iblock)= &
1164 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1165 ggthet(lll,llll,ll,i,j,k,iblock)= &
1166 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1175 ! Control printout of the coefficients of virtual-bond-angle potentials
1178 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1183 write (iout,'(//4a)') &
1184 'Type ',onelett(i),onelett(j),onelett(k)
1185 write (iout,'(//a,10x,a)') " l","a[l]"
1186 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1187 write (iout,'(i2,1pe15.5)') &
1188 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1190 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1191 "b",l,"c",l,"d",l,"e",l
1193 write (iout,'(i2,4(1pe15.5))') m,&
1194 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1195 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1199 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1200 "f+",l,"f-",l,"g+",l,"g-",l
1203 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1204 ffthet(n,m,l,i,j,k,iblock),&
1205 ffthet(m,n,l,i,j,k,iblock),&
1206 ggthet(n,m,l,i,j,k,iblock),&
1207 ggthet(m,n,l,i,j,k,iblock)
1218 !C here will be the apropriate recalibrating for D-aminoacid
1219 read (ithep,*,end=121,err=121) nthetyp
1220 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1221 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1222 do i=-nthetyp+1,nthetyp-1
1223 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1224 do j=0,nbend_kcc_Tb(i)
1225 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1229 write (iout,'(a)') &
1230 "Parameters of the valence-only potentials"
1231 do i=-nthetyp+1,nthetyp-1
1232 write (iout,'(2a)') "Type ",toronelet(i)
1233 do k=0,nbend_kcc_Tb(i)
1234 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1240 write (2,*) "Start reading THETA_PDB",ithep_pdb
1242 ! write (2,*) 'i=',i
1243 read (ithep_pdb,*,err=111,end=111) &
1244 a0thet(i),(athet(j,i,1,1),j=1,2),&
1245 (bthet(j,i,1,1),j=1,2)
1246 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1247 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1248 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1249 sigc0(i)=sigc0(i)**2
1252 athet(1,i,1,-1)=athet(1,i,1,1)
1253 athet(2,i,1,-1)=athet(2,i,1,1)
1254 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1255 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1256 athet(1,i,-1,1)=-athet(1,i,1,1)
1257 athet(2,i,-1,1)=-athet(2,i,1,1)
1258 bthet(1,i,-1,1)=bthet(1,i,1,1)
1259 bthet(2,i,-1,1)=bthet(2,i,1,1)
1262 a0thet(i)=a0thet(-i)
1263 athet(1,i,-1,-1)=athet(1,-i,1,1)
1264 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1265 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1266 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1267 athet(1,i,-1,1)=athet(1,-i,1,1)
1268 athet(2,i,-1,1)=-athet(2,-i,1,1)
1269 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1270 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1271 athet(1,i,1,-1)=-athet(1,-i,1,1)
1272 athet(2,i,1,-1)=athet(2,-i,1,1)
1273 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1274 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1275 theta0(i)=theta0(-i)
1279 polthet(j,i)=polthet(j,-i)
1282 gthet(j,i)=gthet(j,-i)
1285 write (2,*) "End reading THETA_PDB"
1289 !--------------- Reading theta parameters for nucleic acid-------
1290 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1291 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1292 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1293 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1294 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1295 nthetyp_nucl+1,nthetyp_nucl+1))
1296 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1297 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1298 nthetyp_nucl+1,nthetyp_nucl+1))
1299 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1300 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1301 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1302 nthetyp_nucl+1,nthetyp_nucl+1))
1303 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1304 nthetyp_nucl+1,nthetyp_nucl+1))
1305 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1306 nthetyp_nucl+1,nthetyp_nucl+1))
1307 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1308 nthetyp_nucl+1,nthetyp_nucl+1))
1309 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1310 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1311 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1312 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1313 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1314 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1316 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1317 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1319 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1321 aa0thet_nucl(:,:,:)=0.0d0
1322 aathet_nucl(:,:,:,:)=0.0d0
1323 bbthet_nucl(:,:,:,:,:)=0.0d0
1324 ccthet_nucl(:,:,:,:,:)=0.0d0
1325 ddthet_nucl(:,:,:,:,:)=0.0d0
1326 eethet_nucl(:,:,:,:,:)=0.0d0
1327 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1328 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1333 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1334 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1335 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1336 read (ithep_nucl,*,end=111,err=111) &
1337 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1338 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1339 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1340 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1341 read (ithep_nucl,*,end=111,err=111) &
1342 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1343 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1344 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1349 !-------------------------------------------
1350 allocate(nlob(ntyp1)) !(ntyp1)
1351 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1352 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1353 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1360 gaussc(:,:,:,:)=0.0D0
1364 ! Read the parameters of the probability distribution/energy expression
1365 ! of the side chains.
1368 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1372 dsc_inv(i)=1.0D0/dsc(i)
1383 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1384 ((blower(k,l,1),l=1,k),k=1,3)
1385 censc(1,1,-i)=censc(1,1,i)
1386 censc(2,1,-i)=censc(2,1,i)
1387 censc(3,1,-i)=-censc(3,1,i)
1389 read (irotam,*,end=112,err=112) bsc(j,i)
1390 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1391 ((blower(k,l,j),l=1,k),k=1,3)
1392 censc(1,j,-i)=censc(1,j,i)
1393 censc(2,j,-i)=censc(2,j,i)
1394 censc(3,j,-i)=-censc(3,j,i)
1395 ! BSC is amplitude of Gaussian
1402 akl=akl+blower(k,m,j)*blower(l,m,j)
1406 if (((k.eq.3).and.(l.ne.3)) &
1407 .or.((l.eq.3).and.(k.ne.3))) then
1408 gaussc(k,l,j,-i)=-akl
1409 gaussc(l,k,j,-i)=-akl
1411 gaussc(k,l,j,-i)=akl
1412 gaussc(l,k,j,-i)=akl
1421 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1424 if (nlobi.gt.0) then
1426 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1427 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1428 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1429 'log h',(bsc(j,i),j=1,nlobi)
1430 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1431 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1433 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1434 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1437 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1438 write (iout,'(a,f10.4,4(16x,f10.4))') &
1439 'Center ',(bsc(j,i),j=1,nlobi)
1440 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1449 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1450 ! added by Urszula Kozlowska 07/11/2007
1452 !el Maximum number of SC local term fitting function coefficiants
1453 !el integer,parameter :: maxsccoef=65
1455 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1458 read (irotam,*,end=112,err=112)
1460 read (irotam,*,end=112,err=112)
1463 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1467 !---------reading nucleic acid parameters for rotamers-------------------
1468 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1469 do i=1,ntyp_molec(2)
1470 read (irotam_nucl,*,end=112,err=112)
1472 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1478 write (iout,*) "Base rotamer parameters"
1479 do i=1,ntyp_molec(2)
1480 write (iout,'(a)') restyp(i,2)
1481 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1486 ! Read the parameters of the probability distribution/energy expression
1487 ! of the side chains.
1489 write (2,*) "Start reading ROTAM_PDB"
1491 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1495 dsc_inv(i)=1.0D0/dsc(i)
1506 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1507 ((blower(k,l,1),l=1,k),k=1,3)
1509 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1510 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1511 ((blower(k,l,j),l=1,k),k=1,3)
1518 akl=akl+blower(k,m,j)*blower(l,m,j)
1528 write (2,*) "End reading ROTAM_PDB"
1534 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1535 !C interaction energy of the Gly, Ala, and Pro prototypes.
1537 read (ifourier,*) nloctyp
1538 SPLIT_FOURIERTOR = nloctyp.lt.0
1539 nloctyp = iabs(nloctyp)
1540 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1541 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1542 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1543 !C allocate(ctilde(2,2,nres))
1544 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1545 !C allocate(gtb1(2,nres))
1546 !C allocate(gtb2(2,nres))
1547 !C allocate(cc(2,2,nres))
1548 !C allocate(dd(2,2,nres))
1549 !C allocate(ee(2,2,nres))
1550 !C allocate(gtcc(2,2,nres))
1551 !C allocate(gtdd(2,2,nres))
1552 !C allocate(gtee(2,2,nres))
1555 allocate(itype2loc(-ntyp1:ntyp1))
1556 allocate(iloctyp(-nloctyp:nloctyp))
1557 allocate(bnew1(3,2,-nloctyp:nloctyp))
1558 allocate(bnew2(3,2,-nloctyp:nloctyp))
1559 allocate(ccnew(3,2,-nloctyp:nloctyp))
1560 allocate(ddnew(3,2,-nloctyp:nloctyp))
1561 allocate(e0new(3,-nloctyp:nloctyp))
1562 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1563 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1564 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1565 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1566 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1567 allocate(e0newtor(3,-nloctyp:nloctyp))
1568 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1570 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1571 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1572 itype2loc(ntyp1)=nloctyp
1573 iloctyp(nloctyp)=ntyp1
1575 itype2loc(-i)=-itype2loc(i)
1578 allocate(iloctyp(-nloctyp:nloctyp))
1579 allocate(itype2loc(-ntyp1:ntyp1))
1586 iloctyp(-i)=-iloctyp(i)
1588 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1589 !c write (iout,*) "nloctyp",nloctyp,
1590 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1591 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1592 !c write (iout,*) "nloctyp",nloctyp,
1593 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1596 !c write (iout,*) "NEWCORR",i
1597 read (ifourier,*,end=115,err=115)
1600 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1603 !c write (iout,*) "NEWCORR BNEW1"
1604 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1607 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1610 !c write (iout,*) "NEWCORR BNEW2"
1611 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1613 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1614 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1616 !c write (iout,*) "NEWCORR CCNEW"
1617 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1619 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1620 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1622 !c write (iout,*) "NEWCORR DDNEW"
1623 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1627 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1631 !c write (iout,*) "NEWCORR EENEW1"
1632 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1634 read (ifourier,*,end=115,err=115) e0new(ii,i)
1636 !c write (iout,*) (e0new(ii,i),ii=1,3)
1638 !c write (iout,*) "NEWCORR EENEW"
1641 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1642 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1643 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1644 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1649 bnew1(ii,1,-i)= bnew1(ii,1,i)
1650 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1651 bnew2(ii,1,-i)= bnew2(ii,1,i)
1652 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1655 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1656 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1657 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1658 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1659 ccnew(ii,1,-i)=ccnew(ii,1,i)
1660 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1661 ddnew(ii,1,-i)=ddnew(ii,1,i)
1662 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1664 e0new(1,-i)= e0new(1,i)
1665 e0new(2,-i)=-e0new(2,i)
1666 e0new(3,-i)=-e0new(3,i)
1668 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1669 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1670 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1671 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1675 write (iout,'(a)') "Coefficients of the multibody terms"
1676 do i=-nloctyp+1,nloctyp-1
1677 write (iout,*) "Type: ",onelet(iloctyp(i))
1678 write (iout,*) "Coefficients of the expansion of B1"
1680 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1682 write (iout,*) "Coefficients of the expansion of B2"
1684 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1686 write (iout,*) "Coefficients of the expansion of C"
1687 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1688 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1689 write (iout,*) "Coefficients of the expansion of D"
1690 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1691 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1692 write (iout,*) "Coefficients of the expansion of E"
1693 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1696 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1701 IF (SPLIT_FOURIERTOR) THEN
1703 !c write (iout,*) "NEWCORR TOR",i
1704 read (ifourier,*,end=115,err=115)
1707 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1710 !c write (iout,*) "NEWCORR BNEW1 TOR"
1711 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1714 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1717 !c write (iout,*) "NEWCORR BNEW2 TOR"
1718 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1720 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1721 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1723 !c write (iout,*) "NEWCORR CCNEW TOR"
1724 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1726 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1727 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1729 !c write (iout,*) "NEWCORR DDNEW TOR"
1730 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1734 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1738 !c write (iout,*) "NEWCORR EENEW1 TOR"
1739 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1741 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1743 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1745 !c write (iout,*) "NEWCORR EENEW TOR"
1748 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1749 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1750 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1751 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1756 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1757 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1758 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1759 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1762 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1763 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1764 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1765 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1767 e0newtor(1,-i)= e0newtor(1,i)
1768 e0newtor(2,-i)=-e0newtor(2,i)
1769 e0newtor(3,-i)=-e0newtor(3,i)
1771 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1772 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1773 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1774 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1778 write (iout,'(a)') &
1779 "Single-body coefficients of the torsional potentials"
1780 do i=-nloctyp+1,nloctyp-1
1781 write (iout,*) "Type: ",onelet(iloctyp(i))
1782 write (iout,*) "Coefficients of the expansion of B1tor"
1784 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1785 j,(bnew1tor(k,j,i),k=1,3)
1787 write (iout,*) "Coefficients of the expansion of B2tor"
1789 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1790 j,(bnew2tor(k,j,i),k=1,3)
1792 write (iout,*) "Coefficients of the expansion of Ctor"
1793 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1794 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1795 write (iout,*) "Coefficients of the expansion of Dtor"
1796 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1797 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1798 write (iout,*) "Coefficients of the expansion of Etor"
1799 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1802 write (iout,'(1hE,2i1,2f10.5)') &
1803 j,k,(eenewtor(l,j,k,i),l=1,2)
1809 do i=-nloctyp+1,nloctyp-1
1812 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1817 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1821 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1822 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1823 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1824 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1827 ENDIF !SPLIT_FOURIER_TOR
1829 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1830 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1831 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1832 allocate(b(13,-nloctyp-1:nloctyp+1))
1834 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1836 read (ifourier,*,end=115,err=115)
1837 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1839 write (iout,*) 'Type ',onelet(iloctyp(i))
1840 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1854 !c B1tilde(1,i) = b(3)
1855 !c! B1tilde(2,i) =-b(5)
1856 !c! B1tilde(1,-i) =-b(3)
1857 !c! B1tilde(2,-i) =b(5)
1858 !c! b1tilde(1,i)=0.0d0
1859 !c b1tilde(2,i)=0.0d0
1864 !cc B1tilde(1,i) = b(3,i)
1865 !cc B1tilde(2,i) =-b(5,i)
1866 !c B1tilde(1,-i) =-b(3,i)
1867 !c B1tilde(2,-i) =b(5,i)
1868 !cc b1tilde(1,i)=0.0d0
1869 !cc b1tilde(2,i)=0.0d0
1870 !cc B2(1,i) = b(2,i)
1871 !cc B2(2,i) = b(4,i)
1873 !c B2(2,-i) =-b(4,i)
1877 CCold(1,1,i)= b(7,i)
1878 CCold(2,2,i)=-b(7,i)
1879 CCold(2,1,i)= b(9,i)
1880 CCold(1,2,i)= b(9,i)
1881 CCold(1,1,-i)= b(7,i)
1882 CCold(2,2,-i)=-b(7,i)
1883 CCold(2,1,-i)=-b(9,i)
1884 CCold(1,2,-i)=-b(9,i)
1889 !c Ctilde(1,1,i)= CCold(1,1,i)
1890 !c Ctilde(1,2,i)= CCold(1,2,i)
1891 !c Ctilde(2,1,i)=-CCold(2,1,i)
1892 !c Ctilde(2,2,i)=-CCold(2,2,i)
1897 !c Ctilde(1,1,i)= CCold(1,1,i)
1898 !c Ctilde(1,2,i)= CCold(1,2,i)
1899 !c Ctilde(2,1,i)=-CCold(2,1,i)
1900 !c Ctilde(2,2,i)=-CCold(2,2,i)
1902 !c Ctilde(1,1,i)=0.0d0
1903 !c Ctilde(1,2,i)=0.0d0
1904 !c Ctilde(2,1,i)=0.0d0
1905 !c Ctilde(2,2,i)=0.0d0
1906 DDold(1,1,i)= b(6,i)
1907 DDold(2,2,i)=-b(6,i)
1908 DDold(2,1,i)= b(8,i)
1909 DDold(1,2,i)= b(8,i)
1910 DDold(1,1,-i)= b(6,i)
1911 DDold(2,2,-i)=-b(6,i)
1912 DDold(2,1,-i)=-b(8,i)
1913 DDold(1,2,-i)=-b(8,i)
1918 !c Dtilde(1,1,i)= DD(1,1,i)
1919 !c Dtilde(1,2,i)= DD(1,2,i)
1920 !c Dtilde(2,1,i)=-DD(2,1,i)
1921 !c Dtilde(2,2,i)=-DD(2,2,i)
1923 !c Dtilde(1,1,i)=0.0d0
1924 !c Dtilde(1,2,i)=0.0d0
1925 !c Dtilde(2,1,i)=0.0d0
1926 !c Dtilde(2,2,i)=0.0d0
1927 EEold(1,1,i)= b(10,i)+b(11,i)
1928 EEold(2,2,i)=-b(10,i)+b(11,i)
1929 EEold(2,1,i)= b(12,i)-b(13,i)
1930 EEold(1,2,i)= b(12,i)+b(13,i)
1931 EEold(1,1,-i)= b(10,i)+b(11,i)
1932 EEold(2,2,-i)=-b(10,i)+b(11,i)
1933 EEold(2,1,-i)=-b(12,i)+b(13,i)
1934 EEold(1,2,-i)=-b(12,i)-b(13,i)
1935 write(iout,*) "TU DOCHODZE"
1941 !c ee(2,1,i)=ee(1,2,i)
1946 "Coefficients of the cumulants (independent of valence angles)"
1947 do i=-nloctyp+1,nloctyp-1
1948 write (iout,*) 'Type ',onelet(iloctyp(i))
1950 write(iout,'(2f10.5)') B(3,i),B(5,i)
1952 write(iout,'(2f10.5)') B(2,i),B(4,i)
1955 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1959 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1963 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1972 ! Read torsional parameters in old format
1974 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1976 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1977 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1978 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1980 !el from energy module--------
1981 allocate(v1(nterm_old,ntortyp,ntortyp))
1982 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1983 !el---------------------------
1988 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1994 write (iout,'(/a/)') 'Torsional constants:'
1997 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1998 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2004 ! Read torsional parameters
2006 IF (TOR_MODE.eq.0) THEN
2007 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2008 read (itorp,*,end=113,err=113) ntortyp
2009 !el from energy module---------
2010 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2011 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2013 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2014 allocate(vlor2(maxlor,ntortyp,ntortyp))
2015 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2016 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2018 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2019 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2020 !el---------------------------
2023 !el---------------------------
2025 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2027 itortyp(i)=-itortyp(-i)
2029 itortyp(ntyp1)=ntortyp
2030 itortyp(-ntyp1)=-ntortyp
2032 write (iout,*) 'ntortyp',ntortyp
2034 do j=-ntortyp+1,ntortyp-1
2035 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2037 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2038 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2041 do k=1,nterm(i,j,iblock)
2042 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2044 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2045 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2046 v0ij=v0ij+si*v1(k,i,j,iblock)
2048 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2049 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2050 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2052 do k=1,nlor(i,j,iblock)
2053 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2054 vlor2(k,i,j),vlor3(k,i,j)
2055 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2058 v0(-i,-j,iblock)=v0ij
2064 write (iout,'(/a/)') 'Torsional constants:'
2066 do i=-ntortyp,ntortyp
2067 do j=-ntortyp,ntortyp
2068 write (iout,*) 'ityp',i,' jtyp',j
2069 write (iout,*) 'Fourier constants'
2070 do k=1,nterm(i,j,iblock)
2071 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2074 write (iout,*) 'Lorenz constants'
2075 do k=1,nlor(i,j,iblock)
2076 write (iout,'(3(1pe15.5))') &
2077 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2083 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2085 ! 6/23/01 Read parameters for double torsionals
2087 !el from energy module------------
2088 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2089 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2090 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2091 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2092 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2093 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2094 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2095 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2096 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2097 !---------------------------------
2101 do j=-ntortyp+1,ntortyp-1
2102 do k=-ntortyp+1,ntortyp-1
2103 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2104 ! write (iout,*) "OK onelett",
2107 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2108 .or. t3.ne.toronelet(k)) then
2109 write (iout,*) "Error in double torsional parameter file",&
2112 call MPI_Finalize(Ierror)
2114 stop "Error in double torsional parameter file"
2116 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2117 ntermd_2(i,j,k,iblock)
2118 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2119 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2120 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2121 ntermd_1(i,j,k,iblock))
2122 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2123 ntermd_1(i,j,k,iblock))
2124 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2125 ntermd_1(i,j,k,iblock))
2126 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2127 ntermd_1(i,j,k,iblock))
2128 ! Martix of D parameters for one dimesional foureir series
2129 do l=1,ntermd_1(i,j,k,iblock)
2130 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2131 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2132 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2133 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2134 ! write(iout,*) "whcodze" ,
2135 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2137 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2138 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2139 v2s(m,l,i,j,k,iblock),&
2140 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2141 ! Martix of D parameters for two dimesional fourier series
2142 do l=1,ntermd_2(i,j,k,iblock)
2144 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2145 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2146 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2147 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2156 write (iout,*) 'Constants for double torsionals'
2159 do j=-ntortyp+1,ntortyp-1
2160 do k=-ntortyp+1,ntortyp-1
2161 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2162 ' nsingle',ntermd_1(i,j,k,iblock),&
2163 ' ndouble',ntermd_2(i,j,k,iblock)
2165 write (iout,*) 'Single angles:'
2166 do l=1,ntermd_1(i,j,k,iblock)
2167 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2168 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2169 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2170 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2173 write (iout,*) 'Pairs of angles:'
2174 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2175 do l=1,ntermd_2(i,j,k,iblock)
2176 write (iout,'(i5,20f10.5)') &
2177 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2180 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2181 do l=1,ntermd_2(i,j,k,iblock)
2182 write (iout,'(i5,20f10.5)') &
2183 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2184 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2192 ELSE IF (TOR_MODE.eq.1) THEN
2194 !C read valence-torsional parameters
2195 read (itorp,*,end=121,err=121) ntortyp
2197 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2199 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2200 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2203 itype2loc(i)=itortyp(i)
2207 itortyp(i)=-itortyp(-i)
2209 do i=-ntortyp+1,ntortyp-1
2210 do j=-ntortyp+1,ntortyp-1
2211 !C first we read the cos and sin gamma parameters
2212 read (itorp,'(13x,a)',end=121,err=121) string
2213 write (iout,*) i,j,string
2214 read (itorp,*,end=121,err=121) &
2215 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2216 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2217 do k=1,nterm_kcc(j,i)
2218 do l=1,nterm_kcc_Tb(j,i)
2219 do ll=1,nterm_kcc_Tb(j,i)
2220 read (itorp,*,end=121,err=121) ii,jj,kk, &
2221 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2229 !c AL 4/8/16: Calculate coefficients from one-body parameters
2231 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2232 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2233 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2234 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2235 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2238 print *,i,itortyp(i)
2239 itortyp(i)=itype2loc(i)
2242 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2244 do i=-ntortyp+1,ntortyp-1
2245 do j=-ntortyp+1,ntortyp-1
2248 do k=1,nterm_kcc_Tb(j,i)
2249 do l=1,nterm_kcc_Tb(j,i)
2250 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2251 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2252 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2253 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2256 do k=1,nterm_kcc_Tb(j,i)
2257 do l=1,nterm_kcc_Tb(j,i)
2259 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2260 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2261 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2262 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2264 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2265 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2266 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2267 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2271 !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)
2275 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2280 if (tor_mode.gt.0 .and. lprint) then
2281 !c Print valence-torsional parameters
2282 write (iout,'(a)') &
2283 "Parameters of the valence-torsional potentials"
2284 do i=-ntortyp+1,ntortyp-1
2285 do j=-ntortyp+1,ntortyp-1
2286 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2287 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2288 do k=1,nterm_kcc(j,i)
2289 do l=1,nterm_kcc_Tb(j,i)
2290 do ll=1,nterm_kcc_Tb(j,i)
2291 write (iout,'(3i5,2f15.4)')&
2292 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2300 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2301 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2302 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2303 !el from energy module---------
2304 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2305 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2307 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2308 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2309 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2310 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2312 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2313 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2314 !el---------------------------
2317 !el--------------------
2318 read (itorp_nucl,*,end=113,err=113) &
2319 (itortyp_nucl(i),i=1,ntyp_molec(2))
2320 ! print *,itortyp_nucl(:)
2321 !c write (iout,*) 'ntortyp',ntortyp
2324 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2325 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2328 do k=1,nterm_nucl(i,j)
2329 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2330 v0ij=v0ij+si*v1_nucl(k,i,j)
2333 do k=1,nlor_nucl(i,j)
2334 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2335 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2336 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2342 ! Read of Side-chain backbone correlation parameters
2343 ! Modified 11 May 2012 by Adasko
2346 read (isccor,*,end=119,err=119) nsccortyp
2348 !el from module energy-------------
2349 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2350 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2351 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2352 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2353 !-----------------------------------
2355 !el from module energy-------------
2356 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2358 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2360 isccortyp(i)=-isccortyp(-i)
2362 iscprol=isccortyp(20)
2363 ! write (iout,*) 'ntortyp',ntortyp
2365 !c maxinter is maximum interaction sites
2366 !el from module energy---------
2367 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2368 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2369 -nsccortyp:nsccortyp))
2370 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2371 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2372 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2373 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2374 !-----------------------------------
2378 read (isccor,*,end=119,err=119) &
2379 nterm_sccor(i,j),nlor_sccor(i,j)
2385 nterm_sccor(-i,j)=nterm_sccor(i,j)
2386 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2387 nterm_sccor(i,-j)=nterm_sccor(i,j)
2388 do k=1,nterm_sccor(i,j)
2389 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2391 if (j.eq.iscprol) then
2392 if (i.eq.isccortyp(10)) then
2393 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2394 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2396 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2397 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2398 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2399 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2400 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2401 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2402 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2403 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2406 if (i.eq.isccortyp(10)) then
2407 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2408 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2410 if (j.eq.isccortyp(10)) then
2411 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2412 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2414 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2415 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2416 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2417 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2418 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2419 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2423 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2424 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2425 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2426 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2429 do k=1,nlor_sccor(i,j)
2430 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2431 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2432 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2433 (1+vlor3sccor(k,i,j)**2)
2435 v0sccor(l,i,j)=v0ijsccor
2436 v0sccor(l,-i,j)=v0ijsccor1
2437 v0sccor(l,i,-j)=v0ijsccor2
2438 v0sccor(l,-i,-j)=v0ijsccor3
2444 !el from module energy-------------
2445 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2447 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2448 ! write (iout,*) 'ntortyp',ntortyp
2450 !c maxinter is maximum interaction sites
2451 !el from module energy---------
2452 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2453 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2454 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2455 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2456 !-----------------------------------
2460 read (isccor,*,end=119,err=119) &
2461 nterm_sccor(i,j),nlor_sccor(i,j)
2465 do k=1,nterm_sccor(i,j)
2466 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2468 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2471 do k=1,nlor_sccor(i,j)
2472 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2473 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2474 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2475 (1+vlor3sccor(k,i,j)**2)
2477 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2486 write (iout,'(/a/)') 'Torsional constants:'
2489 write (iout,*) 'ityp',i,' jtyp',j
2490 write (iout,*) 'Fourier constants'
2491 do k=1,nterm_sccor(i,j)
2492 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2494 write (iout,*) 'Lorenz constants'
2495 do k=1,nlor_sccor(i,j)
2496 write (iout,'(3(1pe15.5))') &
2497 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2504 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2505 ! interaction energy of the Gly, Ala, and Pro prototypes.
2508 ! Read electrostatic-interaction parameters
2513 write (iout,'(/a)') 'Electrostatic interaction constants:'
2514 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2515 'IT','JT','APP','BPP','AEL6','AEL3'
2517 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2518 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2519 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2520 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2525 app (i,j)=epp(i,j)*rri*rri
2526 bpp (i,j)=-2.0D0*epp(i,j)*rri
2527 ael6(i,j)=elpp6(i,j)*4.2D0**6
2528 ael3(i,j)=elpp3(i,j)*4.2D0**3
2530 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2536 ! Read side-chain interaction parameters.
2538 !el from module energy - COMMON.INTERACT-------
2539 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2540 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2541 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2543 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2544 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2545 allocate(epslip(ntyp,ntyp))
2553 !--------------------------------
2555 read (isidep,*,end=117,err=117) ipot,expon
2556 if (ipot.lt.1 .or. ipot.gt.5) then
2557 write (iout,'(2a)') 'Error while reading SC interaction',&
2558 'potential file - unknown potential type.'
2560 call MPI_Finalize(Ierror)
2566 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2567 ', exponents are ',expon,2*expon
2568 ! goto (10,20,30,30,40) ipot
2570 !----------------------- LJ potential ---------------------------------
2572 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2573 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2574 (sigma0(i),i=1,ntyp)
2576 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2577 write (iout,'(a/)') 'The epsilon array:'
2578 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2579 write (iout,'(/a)') 'One-body parameters:'
2580 write (iout,'(a,4x,a)') 'residue','sigma'
2581 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2584 !----------------------- LJK potential --------------------------------
2586 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2587 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2588 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2590 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2591 write (iout,'(a/)') 'The epsilon array:'
2592 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2593 write (iout,'(/a)') 'One-body parameters:'
2594 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2595 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2599 !---------------------- GB or BP potential -----------------------------
2602 ! print *,"I AM in SCELE",scelemode
2603 if (scelemode.eq.0) then
2605 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2607 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2608 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2609 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2610 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2612 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2615 ! For the GB potential convert sigma'**2 into chi'
2618 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2622 write (iout,'(/a/)') 'Parameters of the BP potential:'
2623 write (iout,'(a/)') 'The epsilon array:'
2624 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2625 write (iout,'(/a)') 'One-body parameters:'
2626 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2628 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2629 chip(i),alp(i),i=1,ntyp)
2632 ! print *,ntyp,"NTYP"
2633 allocate(icharge(ntyp1))
2634 ! print *,ntyp,icharge(i)
2636 read (isidep,*) (icharge(i),i=1,ntyp)
2637 print *,ntyp,icharge(i)
2638 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2639 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2640 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2641 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2642 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2643 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2644 allocate(epsintab(ntyp,ntyp))
2645 allocate(dtail(2,ntyp,ntyp))
2646 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2647 allocate(wqdip(2,ntyp,ntyp))
2648 allocate(wstate(4,ntyp,ntyp))
2649 allocate(dhead(2,2,ntyp,ntyp))
2650 allocate(nstate(ntyp,ntyp))
2651 allocate(debaykap(ntyp,ntyp))
2653 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2654 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2658 ! write (*,*) "Im in ALAB", i, " ", j
2660 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2661 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2662 chis(i,j),chis(j,i), &
2663 nstate(i,j),(wstate(k,i,j),k=1,4), &
2664 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2665 dtail(1,i,j),dtail(2,i,j), &
2666 epshead(i,j),sig0head(i,j), &
2667 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2668 alphapol(i,j),alphapol(j,i), &
2669 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
2670 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2676 sigma(i,j) = sigma(j,i)
2677 sigmap1(i,j)=sigmap1(j,i)
2678 sigmap2(i,j)=sigmap2(j,i)
2679 sigiso1(i,j)=sigiso1(j,i)
2680 sigiso2(i,j)=sigiso2(j,i)
2681 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2682 nstate(i,j) = nstate(j,i)
2683 dtail(1,i,j) = dtail(1,j,i)
2684 dtail(2,i,j) = dtail(2,j,i)
2686 alphasur(k,i,j) = alphasur(k,j,i)
2687 wstate(k,i,j) = wstate(k,j,i)
2688 alphiso(k,i,j) = alphiso(k,j,i)
2691 dhead(2,1,i,j) = dhead(1,1,j,i)
2692 dhead(2,2,i,j) = dhead(1,2,j,i)
2693 dhead(1,1,i,j) = dhead(2,1,j,i)
2694 dhead(1,2,i,j) = dhead(2,2,j,i)
2696 epshead(i,j) = epshead(j,i)
2697 sig0head(i,j) = sig0head(j,i)
2700 wqdip(k,i,j) = wqdip(k,j,i)
2703 wquad(i,j) = wquad(j,i)
2704 epsintab(i,j) = epsintab(j,i)
2705 debaykap(i,j)=debaykap(j,i)
2706 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2713 !--------------------- GBV potential -----------------------------------
2715 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2716 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2717 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2718 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2720 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2721 write (iout,'(a/)') 'The epsilon array:'
2722 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2723 write (iout,'(/a)') 'One-body parameters:'
2724 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2725 's||/s_|_^2',' chip ',' alph '
2726 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2727 sigii(i),chip(i),alp(i),i=1,ntyp)
2730 write(iout,*)"Wrong ipot"
2736 !-----------------------------------------------------------------------
2737 ! Calculate the "working" parameters of SC interactions.
2739 !el from module energy - COMMON.INTERACT-------
2740 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2741 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2742 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2743 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2744 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2745 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2747 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2753 if (scelemode.eq.0) then
2762 sc_aa_tube_par(:)=0.0d0
2763 sc_bb_tube_par(:)=0.0d0
2765 !--------------------------------
2770 epslip(i,j)=epslip(j,i)
2773 if (scelemode.eq.0) then
2776 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2777 sigma(j,i)=sigma(i,j)
2778 rs0(i,j)=dwa16*sigma(i,j)
2783 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2784 'Working parameters of the SC interactions:',&
2785 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2790 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2792 ! print *,"SIGMA ZLA?",sigma(i,j)
2800 sigeps=dsign(1.0D0,epsij)
2802 aa_aq(i,j)=epsij*rrij*rrij
2803 ! print *,"ADASKO",epsij,rrij,expon
2804 bb_aq(i,j)=-sigeps*epsij*rrij
2805 aa_aq(j,i)=aa_aq(i,j)
2806 bb_aq(j,i)=bb_aq(i,j)
2807 epsijlip=epslip(i,j)
2808 sigeps=dsign(1.0D0,epsijlip)
2809 epsijlip=dabs(epsijlip)
2810 aa_lip(i,j)=epsijlip*rrij*rrij
2811 bb_lip(i,j)=-sigeps*epsijlip*rrij
2812 aa_lip(j,i)=aa_lip(i,j)
2813 bb_lip(j,i)=bb_lip(i,j)
2814 !C write(iout,*) aa_lip
2815 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2816 sigt1sq=sigma0(i)**2
2817 sigt2sq=sigma0(j)**2
2820 ratsig1=sigt2sq/sigt1sq
2821 ratsig2=1.0D0/ratsig1
2822 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2823 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2824 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2828 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2829 sigmaii(i,j)=rsum_max
2830 sigmaii(j,i)=rsum_max
2832 ! sigmaii(i,j)=r0(i,j)
2833 ! sigmaii(j,i)=r0(i,j)
2835 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2836 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2837 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2838 augm(i,j)=epsij*r_augm**(2*expon)
2839 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2846 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2847 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2848 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2853 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2854 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2855 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2856 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2857 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2858 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2859 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2860 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2861 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2862 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2863 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2864 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2865 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2866 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2867 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2868 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2877 read (isidep_nucl,*) ipot_nucl
2878 ! print *,"TU?!",ipot_nucl
2879 if (ipot_nucl.eq.1) then
2880 do i=1,ntyp_molec(2)
2881 do j=i,ntyp_molec(2)
2882 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2883 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2887 do i=1,ntyp_molec(2)
2888 do j=i,ntyp_molec(2)
2889 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2890 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2891 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2895 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2896 do i=1,ntyp_molec(2)
2897 do j=i,ntyp_molec(2)
2898 rrij=sigma_nucl(i,j)
2902 epsij=4*eps_nucl(i,j)
2903 sigeps=dsign(1.0D0,epsij)
2905 aa_nucl(i,j)=epsij*rrij*rrij
2906 bb_nucl(i,j)=-sigeps*epsij*rrij
2907 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2908 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2909 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2910 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2911 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2912 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2915 aa_nucl(i,j)=aa_nucl(j,i)
2916 bb_nucl(i,j)=bb_nucl(j,i)
2917 ael3_nucl(i,j)=ael3_nucl(j,i)
2918 ael6_nucl(i,j)=ael6_nucl(j,i)
2919 ael63_nucl(i,j)=ael63_nucl(j,i)
2920 ael32_nucl(i,j)=ael32_nucl(j,i)
2921 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2922 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2923 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2924 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2925 eps_nucl(i,j)=eps_nucl(j,i)
2926 sigma_nucl(i,j)=sigma_nucl(j,i)
2927 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2931 write(iout,*) "tube param"
2932 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2933 ccavtubpep,dcavtubpep,tubetranenepep
2934 sigmapeptube=sigmapeptube**6
2935 sigeps=dsign(1.0D0,epspeptube)
2936 epspeptube=dabs(epspeptube)
2937 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2938 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2939 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2941 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2942 ccavtub(i),dcavtub(i),tubetranene(i)
2943 sigmasctube=sigmasctube**6
2944 sigeps=dsign(1.0D0,epssctube)
2945 epssctube=dabs(epssctube)
2946 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2947 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2948 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2950 !-----------------READING SC BASE POTENTIALS-----------------------------
2951 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2952 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2953 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2954 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2955 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2956 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2957 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2958 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2959 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2960 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2961 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2962 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2963 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2964 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2965 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2966 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2969 do i=1,ntyp_molec(1)
2970 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2971 write (*,*) "Im in ", i, " ", j
2972 read(isidep_scbase,*) &
2973 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2974 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2975 write(*,*) "eps",eps_scbase(i,j)
2976 read(isidep_scbase,*) &
2977 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2978 chis_scbase(i,j,1),chis_scbase(i,j,2)
2979 read(isidep_scbase,*) &
2980 dhead_scbasei(i,j), &
2981 dhead_scbasej(i,j), &
2982 rborn_scbasei(i,j),rborn_scbasej(i,j)
2983 read(isidep_scbase,*) &
2984 (wdipdip_scbase(k,i,j),k=1,3), &
2985 (wqdip_scbase(k,i,j),k=1,2)
2986 read(isidep_scbase,*) &
2987 alphapol_scbase(i,j), &
2988 epsintab_scbase(i,j)
2991 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2992 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2994 do i=1,ntyp_molec(1)
2995 do j=1,ntyp_molec(2)-1
2996 epsij=eps_scbase(i,j)
2997 rrij=sigma_scbase(i,j)
3002 sigeps=dsign(1.0D0,epsij)
3004 aa_scbase(i,j)=epsij*rrij*rrij
3005 bb_scbase(i,j)=-sigeps*epsij*rrij
3008 !-----------------READING PEP BASE POTENTIALS-------------------
3009 allocate(eps_pepbase(ntyp_molec(2)))
3010 allocate(sigma_pepbase(ntyp_molec(2)))
3011 allocate(chi_pepbase(ntyp_molec(2),2))
3012 allocate(chipp_pepbase(ntyp_molec(2),2))
3013 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3014 allocate(sigmap1_pepbase(ntyp_molec(2)))
3015 allocate(sigmap2_pepbase(ntyp_molec(2)))
3016 allocate(chis_pepbase(ntyp_molec(2),2))
3017 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3020 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3021 write (*,*) "Im in ", i, " ", j
3022 read(isidep_pepbase,*) &
3023 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3024 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3025 write(*,*) "eps",eps_pepbase(j)
3026 read(isidep_pepbase,*) &
3027 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3028 chis_pepbase(j,1),chis_pepbase(j,2)
3029 read(isidep_pepbase,*) &
3030 (wdipdip_pepbase(k,j),k=1,3)
3032 allocate(aa_pepbase(ntyp_molec(2)))
3033 allocate(bb_pepbase(ntyp_molec(2)))
3035 do j=1,ntyp_molec(2)-1
3036 epsij=eps_pepbase(j)
3037 rrij=sigma_pepbase(j)
3042 sigeps=dsign(1.0D0,epsij)
3044 aa_pepbase(j)=epsij*rrij*rrij
3045 bb_pepbase(j)=-sigeps*epsij*rrij
3047 !--------------READING SC PHOSPHATE-------------------------------------
3048 allocate(eps_scpho(ntyp_molec(1)))
3049 allocate(sigma_scpho(ntyp_molec(1)))
3050 allocate(chi_scpho(ntyp_molec(1),2))
3051 allocate(chipp_scpho(ntyp_molec(1),2))
3052 allocate(alphasur_scpho(4,ntyp_molec(1)))
3053 allocate(sigmap1_scpho(ntyp_molec(1)))
3054 allocate(sigmap2_scpho(ntyp_molec(1)))
3055 allocate(chis_scpho(ntyp_molec(1),2))
3056 allocate(wqq_scpho(ntyp_molec(1)))
3057 allocate(wqdip_scpho(2,ntyp_molec(1)))
3058 allocate(alphapol_scpho(ntyp_molec(1)))
3059 allocate(epsintab_scpho(ntyp_molec(1)))
3060 allocate(dhead_scphoi(ntyp_molec(1)))
3061 allocate(rborn_scphoi(ntyp_molec(1)))
3062 allocate(rborn_scphoj(ntyp_molec(1)))
3063 allocate(alphi_scpho(ntyp_molec(1)))
3067 do j=1,ntyp_molec(1) ! without U then we will take T for U
3068 write (*,*) "Im in scpho ", i, " ", j
3069 read(isidep_scpho,*) &
3070 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3071 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3072 write(*,*) "eps",eps_scpho(j)
3073 read(isidep_scpho,*) &
3074 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3075 chis_scpho(j,1),chis_scpho(j,2)
3076 read(isidep_scpho,*) &
3077 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3078 read(isidep_scpho,*) &
3079 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3083 allocate(aa_scpho(ntyp_molec(1)))
3084 allocate(bb_scpho(ntyp_molec(1)))
3086 do j=1,ntyp_molec(1)
3093 sigeps=dsign(1.0D0,epsij)
3095 aa_scpho(j)=epsij*rrij*rrij
3096 bb_scpho(j)=-sigeps*epsij*rrij
3100 read(isidep_peppho,*) &
3101 eps_peppho,sigma_peppho
3102 read(isidep_peppho,*) &
3103 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3104 read(isidep_peppho,*) &
3105 (wqdip_peppho(k),k=1,2)
3113 sigeps=dsign(1.0D0,epsij)
3115 aa_peppho=epsij*rrij*rrij
3116 bb_peppho=-sigeps*epsij*rrij
3119 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3124 ! Define the SC-p interaction constants (hard-coded; old style)
3127 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3129 ! aad(i,1)=0.3D0*4.0D0**12
3130 ! Following line for constants currently implemented
3131 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3132 aad(i,1)=1.5D0*4.0D0**12
3133 ! aad(i,1)=0.17D0*5.6D0**12
3135 ! "Soft" SC-p repulsion
3137 ! Following line for constants currently implemented
3138 ! aad(i,1)=0.3D0*4.0D0**6
3139 ! "Hard" SC-p repulsion
3140 bad(i,1)=3.0D0*4.0D0**6
3141 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3150 ! 8/9/01 Read the SC-p interaction constants from file
3153 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3156 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3157 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3158 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3159 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3163 write (iout,*) "Parameters of SC-p interactions:"
3165 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3166 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3171 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3173 do i=1,ntyp_molec(2)
3174 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3176 do i=1,ntyp_molec(2)
3177 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3178 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3180 r0pp=1.12246204830937298142*5.16158
3186 ! Define the constants of the disulfide bridge
3190 ! Old arbitrary potential - commented out.
3195 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3196 ! energy surface of diethyl disulfide.
3197 ! A. Liwo and U. Kozlowska, 11/24/03
3215 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3216 allocate(sigiso1cat(ntyp,ntyp),rborncat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3217 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3218 allocate(chiscat(ntyp,ntyp),wquadcat(ntyp,ntyp),chippcat(ntyp,ntyp))
3219 allocate(epsintabcat(ntyp,ntyp))
3220 allocate(dtailcat(2,ntyp,ntyp))
3221 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3222 allocate(wqdipcat(2,ntyp,ntyp))
3223 allocate(wstatecat(4,ntyp,ntyp))
3224 allocate(dheadcat(2,2,ntyp,ntyp))
3225 allocate(nstatecat(ntyp,ntyp))
3226 allocate(debaykapcat(ntyp,ntyp))
3229 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3230 if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3232 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3233 if (oldion.eq.0) then
3234 do i=1,ntyp_molec(5)
3235 read(iion,*) msc(i,5),restok(i,5)
3236 print *,msc(i,5),restok(i,5)
3241 do j=1,ntyp_molec(5)
3242 ! write (*,*) "Im in ALAB", i, " ", j
3244 epscat(i,j),sigmacat(i,j),chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3245 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3246 chiscat(i,j),chiscat(j,i), &
3247 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3248 dtailcat(1,i,j),dtailcat(2,i,j), &
3249 epsheadcat(i,j),sig0headcat(i,j), &
3251 rborncat(i,j),rborncat(j,i),(wqdipcat(k,i,j),k=1,2), &
3252 alphapolcat(i,j),alphapolcat(j,i), &
3253 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3254 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3261 write (iout,'(/a)') "Disulfide bridge parameters:"
3262 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3263 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3264 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3265 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3268 if (shield_mode.gt.0) then
3269 pi=4.0D0*datan(1.0D0)
3270 !C VSolvSphere the volume of solving sphere
3272 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3273 !C there will be no distinction between proline peptide group and normal peptide
3274 !C group in case of shielding parameters
3275 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3276 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3277 write (iout,*) VSolvSphere,VSolvSphere_div
3278 !C long axis of side chain
3280 long_r_sidechain(i)=vbldsc0(1,i)
3281 ! if (scelemode.eq.0) then
3282 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3283 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3285 ! short_r_sidechain(i)=sigma(i,i)
3287 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3294 111 write (iout,*) "Error reading bending energy parameters."
3296 112 write (iout,*) "Error reading rotamer energy parameters."
3298 113 write (iout,*) "Error reading torsional energy parameters."
3300 114 write (iout,*) "Error reading double torsional energy parameters."
3302 115 write (iout,*) &
3303 "Error reading cumulant (multibody energy) parameters."
3305 116 write (iout,*) "Error reading electrostatic energy parameters."
3307 117 write (iout,*) "Error reading side chain interaction parameters."
3309 118 write (iout,*) "Error reading SCp interaction parameters."
3311 119 write (iout,*) "Error reading SCCOR parameters"
3313 121 write (iout,*) "Error in Czybyshev parameters"
3316 call MPI_Finalize(Ierror)
3320 end subroutine parmread
3322 !-----------------------------------------------------------------------------
3324 !-----------------------------------------------------------------------------
3325 subroutine printmat(ldim,m,n,iout,key,a)
3328 character(len=3),dimension(n) :: key
3329 real(kind=8),dimension(ldim,n) :: a
3331 integer :: i,j,k,m,iout,nlim
3335 write (iout,1000) (key(k),k=i,nlim)
3337 1000 format (/5x,8(6x,a3))
3338 1020 format (/80(1h-)/)
3340 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3343 1010 format (a3,2x,8(f9.4))
3345 end subroutine printmat
3346 !-----------------------------------------------------------------------------
3348 !-----------------------------------------------------------------------------
3350 ! Read the PDB file and convert the peptide geometry into virtual-chain
3353 use energy_data, only: itype,istype
3357 ! use control, only: rescode,sugarcode
3358 ! implicit real*8 (a-h,o-z)
3359 ! include 'DIMENSIONS'
3360 ! include 'COMMON.LOCAL'
3361 ! include 'COMMON.VAR'
3362 ! include 'COMMON.CHAIN'
3363 ! include 'COMMON.INTERACT'
3364 ! include 'COMMON.IOUNITS'
3365 ! include 'COMMON.GEO'
3366 ! include 'COMMON.NAMES'
3367 ! include 'COMMON.CONTROL'
3368 ! include 'COMMON.DISTFIT'
3369 ! include 'COMMON.SETUP'
3370 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3372 logical :: lprn=.true.,fail
3373 real(kind=8),dimension(3) :: e1,e2,e3
3374 real(kind=8) :: dcj,efree_temp
3375 character(len=3) :: seq,res,res2
3376 character(len=5) :: atom
3377 character(len=80) :: card
3378 real(kind=8),dimension(3,20) :: sccor
3379 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3380 integer :: isugar,molecprev,firstion
3381 character*1 :: sugar
3383 real(kind=8),dimension(3) :: ccc
3385 integer,dimension(2,maxres/3) :: hfrag_alloc
3386 integer,dimension(4,maxres/3) :: bfrag_alloc
3387 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3388 real(kind=8),dimension(:,:), allocatable :: c_temporary
3389 integer,dimension(:,:) , allocatable :: itype_temporary
3390 integer,dimension(:),allocatable :: istype_temp
3397 ! write (2,*) "UNRES_PDB",unres_pdb
3417 !-----------------------------
3418 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3419 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3420 if(.not. allocated(istype)) allocate(istype(maxres))
3422 read (ipdbin,'(a80)',end=10) card
3423 write (iout,'(a)') card
3424 if (card(:5).eq.'HELIX') then
3427 read(card(22:25),*) hfrag(1,nhfrag)
3428 read(card(34:37),*) hfrag(2,nhfrag)
3430 if (card(:5).eq.'SHEET') then
3433 read(card(24:26),*) bfrag(1,nbfrag)
3434 read(card(35:37),*) bfrag(2,nbfrag)
3435 !rc----------------------------------------
3436 !rc to be corrected !!!
3437 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3438 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3439 !rc----------------------------------------
3441 if (card(:3).eq.'END') then
3443 else if (card(:3).eq.'TER') then
3448 itype(ires_old,molecule)=ntyp1_molec(molecule)
3449 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3450 nres_molec(molecule)=nres_molec(molecule)+2
3452 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3455 dc(j,ires)=sccor(j,iii)
3458 call sccenter(ires,iii,sccor)
3464 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3465 ! Fish out the ATOM cards.
3466 ! write(iout,*) 'card',card(1:20)
3467 ! print *,"ATU ",card(1:6), CARD(3:6)
3469 if (index(card(1:4),'ATOM').gt.0) then
3470 read (card(12:16),*) atom
3471 ! write (iout,*) "! ",atom," !",ires
3472 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3473 read (card(23:26),*) ires
3474 read (card(18:20),'(a3)') res
3475 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3476 ! & " ires_old",ires_old
3477 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3478 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3479 if (ires-ishift+ishift1.ne.ires_old) then
3480 ! Calculate the CM of the preceding residue.
3481 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3483 ! write (iout,*) "Calculating sidechain center iii",iii
3486 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3489 call sccenter(ires_old,iii,sccor)
3493 ! Start new residue.
3494 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3497 else if (ibeg.eq.1) then
3498 write (iout,*) "BEG ires",ires
3500 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3503 nres_molec(molecule)=nres_molec(molecule)+1
3505 ires=ires-ishift+ishift1
3507 ! write (iout,*) "ishift",ishift," ires",ires,&
3508 ! " ires_old",ires_old
3510 else if (ibeg.eq.2) then
3512 ishift=-ires_old+ires-1 !!!!!
3513 ishift1=ishift1-1 !!!!!
3514 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3515 ires=ires-ishift+ishift1
3516 ! print *,ires,ishift,ishift1
3520 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3521 ires=ires-ishift+ishift1
3524 ! print *,'atom',ires,atom
3525 if (res.eq.'ACE' .or. res.eq.'NHE') then
3528 if (atom.eq.'CA '.or.atom.eq.'N ') then
3530 itype(ires,molecule)=rescode(ires,res,0,molecule)
3532 ! nres_molec(molecule)=nres_molec(molecule)+1
3536 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3537 ! nres_molec(molecule)=nres_molec(molecule)+1
3538 read (card(19:19),'(a1)') sugar
3539 isugar=sugarcode(sugar,ires)
3540 ! if (ibeg.eq.1) then
3544 ! print *,"ires=",ires,istype(ires)
3550 ires=ires-ishift+ishift1
3552 ! write (iout,*) "ires_old",ires_old," ires",ires
3553 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3556 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3557 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3558 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3559 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3560 ! print *,ires,ishift,ishift1
3561 ! write (iout,*) "backbone ",atom
3563 write (iout,'(2i3,2x,a,3f8.3)') &
3564 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3567 nres_molec(molecule)=nres_molec(molecule)+1
3569 sccor(j,iii)=c(j,ires)
3571 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3572 atom.eq."C2'" .or. atom.eq."C3'" &
3573 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3574 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3575 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3576 ! print *,ires,ishift,ishift1
3580 ! sccor(j,iii)=c(j,ires)
3583 c(j,ires)=c(j,ires)+ccc(j)/5.0
3585 print *,counter,molecule
3586 if (counter.eq.5) then
3588 nres_molec(molecule)=nres_molec(molecule)+1
3591 ! sccor(j,iii)=c(j,ires)
3595 ! print *, "ATOM",atom(1:3)
3596 else if (atom.eq."C5'") then
3597 read (card(19:19),'(a1)') sugar
3598 isugar=sugarcode(sugar,ires)
3603 ! print *,ires,istype(ires)
3607 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3608 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3609 nres_molec(molecule)=nres_molec(molecule)+1
3610 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3614 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3616 else if ((atom.eq."C1'").and.unres_pdb) then
3618 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3619 ! write (*,*) card(23:27),ires,itype(ires,1)
3620 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3621 atom.ne.'N' .and. atom.ne.'C' .and. &
3622 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3623 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3624 .and. atom.ne.'P '.and. &
3625 atom(1:1).ne.'H' .and. &
3626 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3628 ! write (iout,*) "sidechain ",atom
3629 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3630 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3631 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3633 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3636 ! print *,"IONS",ions,card(1:6)
3637 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3638 if (firstion.eq.0) then
3642 dc(j,ires)=sccor(j,iii)
3645 call sccenter(ires,iii,sccor)
3648 read (card(12:16),*) atom
3649 ! print *,"HETATOM", atom
3650 read (card(18:20),'(a3)') res
3651 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3652 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3653 .or.(atom(1:2).eq.'K ')) &
3656 if (molecule.ne.5) molecprev=molecule
3658 nres_molec(molecule)=nres_molec(molecule)+1
3659 print *,"HERE",nres_molec(molecule)
3661 itype(ires,molecule)=rescode(ires,res,0,molecule)
3662 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3666 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3667 if (ires.eq.0) return
3668 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3671 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3673 nres_molec(molecule)=nres_molec(molecule)-2
3674 print *,'I have',nres, nres_molec(:)
3676 do k=1,4 ! ions are without dummy
3677 if (nres_molec(k).eq.0) cycle
3679 ! write (iout,*) i,itype(i,1)
3680 ! if (itype(i,1).eq.ntyp1) then
3681 ! write (iout,*) "dummy",i,itype(i,1)
3683 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3684 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3688 if (itype(i,k).eq.ntyp1_molec(k)) then
3689 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3690 if (itype(i+2,k).eq.0) then
3691 ! print *,"masz sieczke"
3693 if (itype(i+2,j).ne.0) then
3695 itype(i+1,j)=ntyp1_molec(j)
3696 nres_molec(k)=nres_molec(k)-1
3697 nres_molec(j)=nres_molec(j)+1
3703 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3704 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3705 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3706 ! if (unres_pdb) then
3707 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3708 ! print *,i,'tu dochodze'
3709 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3717 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3721 dcj=(c(j,i-2)-c(j,i-3))/2.0
3722 if (dcj.eq.0) dcj=1.23591524223
3727 else !itype(i+1,1).eq.ntyp1
3728 ! if (unres_pdb) then
3729 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3730 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3737 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3738 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3742 dcj=(c(j,i+3)-c(j,i+2))/2.0
3743 if (dcj.eq.0) dcj=1.23591524223
3748 endif !itype(i+1,1).eq.ntyp1
3749 endif !itype.eq.ntyp1
3753 ! Calculate the CM of the last side chain.
3757 dc(j,ires)=sccor(j,iii)
3760 call sccenter(ires,iii,sccor)
3766 ! print *,"molecule",molecule
3767 if ((itype(nres,1).ne.10)) then
3769 if (molecule.eq.5) molecule=molecprev
3770 itype(nres,molecule)=ntyp1_molec(molecule)
3771 nres_molec(molecule)=nres_molec(molecule)+1
3772 ! if (unres_pdb) then
3773 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3774 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3781 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3785 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3786 c(j,nres)=c(j,nres-1)+dcj
3787 c(j,2*nres)=c(j,nres)
3791 ! print *,'I have',nres, nres_molec(:)
3793 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3795 if (nres.ne.nres0) then
3796 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3798 stop "Error nres value in WHAM input"
3801 !---------------------------------
3802 !el reallocate tables
3805 ! hfrag_alloc(j,i)=hfrag(j,i)
3808 ! bfrag_alloc(j,i)=bfrag(j,i)
3814 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3815 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3816 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3817 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3821 ! hfrag(j,i)=hfrag_alloc(j,i)
3826 ! bfrag(j,i)=bfrag_alloc(j,i)
3829 !el end reallocate tables
3830 !---------------------------------
3838 c(j,2*nres)=c(j,nres)
3841 if (itype(1,1).eq.ntyp1) then
3845 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3846 call refsys(2,3,4,e1,e2,e3,fail)
3853 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3854 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
3858 dcj=(c(j,4)-c(j,3))/2.0
3864 ! First lets assign correct dummy to correct type of chain
3866 if (itype(1,1).eq.ntyp1) then
3867 if (itype(2,1).eq.0) then
3869 if (itype(2,j).ne.0) then
3871 itype(1,j)=ntyp1_molec(j)
3872 nres_molec(1)=nres_molec(1)-1
3873 nres_molec(j)=nres_molec(j)+1
3880 print *,'I have',nres, nres_molec(:)
3882 ! Copy the coordinates to reference coordinates
3888 ! Calculate internal coordinates.
3890 write (iout,'(/a)') &
3891 "Cartesian coordinates of the reference structure"
3892 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3893 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3895 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3896 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3897 (c(j,ires+nres),j=1,3)
3900 ! znamy już nres wiec mozna alokowac tablice
3901 ! Calculate internal coordinates.
3902 if(me.eq.king.or..not.out1file)then
3903 write (iout,'(a)') &
3904 "Backbone and SC coordinates as read from the PDB"
3906 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
3907 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3908 (c(j,nres+ires),j=1,3)
3911 ! NOW LETS ROCK! SORTING
3912 allocate(c_temporary(3,2*nres))
3913 allocate(itype_temporary(nres,5))
3914 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3915 if (.not.allocated(istype)) write(iout,*) &
3916 "SOMETHING WRONG WITH ISTYTPE"
3917 allocate(istype_temp(nres))
3918 itype_temporary(:,:)=0
3922 if (itype(i,k).ne.0) then
3924 c_temporary(j,seqalingbegin)=c(j,i)
3925 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3928 itype_temporary(seqalingbegin,k)=itype(i,k)
3929 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3930 istype_temp(seqalingbegin)=istype(i)
3931 molnum(seqalingbegin)=k
3932 seqalingbegin=seqalingbegin+1
3938 c(j,i)=c_temporary(j,i)
3943 itype(i,k)=itype_temporary(i,k)
3944 istype(i)=istype_temp(i)
3947 ! if (itype(1,1).eq.ntyp1) then
3950 ! if (unres_pdb) then
3951 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3952 ! call refsys(2,3,4,e1,e2,e3,fail)
3959 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3963 ! dcj=(c(j,4)-c(j,3))/2.0
3965 ! c(j,nres+1)=c(j,1)
3971 write (iout,'(/a)') &
3972 "Cartesian coordinates of the reference structure after sorting"
3973 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3974 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3976 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3977 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3978 (c(j,ires+nres),j=1,3)
3982 ! print *,seqalingbegin,nres
3983 if(.not.allocated(vbld)) then
3984 allocate(vbld(2*nres))
3989 if(.not.allocated(vbld_inv)) then
3990 allocate(vbld_inv(2*nres))
3996 if(.not.allocated(theta)) then
3997 allocate(theta(nres+2))
4001 if(.not.allocated(phi)) allocate(phi(nres+2))
4002 if(.not.allocated(alph)) allocate(alph(nres+2))
4003 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4004 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4005 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4006 if(.not.allocated(costtab)) allocate(costtab(nres))
4007 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4008 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4009 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4010 if(.not.allocated(xxref)) allocate(xxref(nres))
4011 if(.not.allocated(yyref)) allocate(yyref(nres))
4012 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4013 if(.not.allocated(dc_norm)) then
4014 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4015 allocate(dc_norm(3,0:2*nres+2))
4019 call int_from_cart(.true.,.false.)
4020 call sc_loc_geom(.false.)
4022 thetaref(i)=theta(i)
4032 dc(j,i)=c(j,i+1)-c(j,i)
4033 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4038 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4039 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4041 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4045 ! Copy the coordinates to reference coordinates
4046 ! Splits to single chain if occurs
4050 ! cref(j,i,cou)=c(j,i)
4054 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4055 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4056 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4057 !-----------------------------
4061 write (iout,*) "symetr", symetr
4064 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4066 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4069 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4074 cref(j,i,cou)=c(j,i)
4075 cref(j,i+nres,cou)=c(j,i+nres)
4077 chain_rep(j,lll,kkk)=c(j,i)
4078 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4082 write (iout,*) chain_length
4083 if (chain_length.eq.0) chain_length=nres
4085 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4086 chain_rep(j,chain_length+nres,symetr) &
4087 =chain_rep(j,chain_length+nres,1)
4090 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4092 ! do kkk=1,chain_length
4093 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4097 ! makes copy of chains
4098 write (iout,*) "symetr", symetr
4103 if (symetr.gt.1) then
4110 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4116 write (iout,*) i,icha
4117 do lll=1,chain_length
4119 if (cou.le.nres) then
4121 kupa=mod(lll,chain_length)
4122 iprzes=(kkk-1)*chain_length+lll
4123 if (kupa.eq.0) kupa=chain_length
4124 write (iout,*) "kupa", kupa
4125 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4126 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4133 !-koniec robienia kopii
4136 write (iout,*) "nowa struktura", nperm
4138 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4140 cref(3,i,kkk),cref(1,nres+i,kkk),&
4141 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4143 100 format (//' alpha-carbon coordinates ',&
4144 ' centroid coordinates'/ &
4145 ' ', 6X,'X',11X,'Y',11X,'Z', &
4146 10X,'X',11X,'Y',11X,'Z')
4147 110 format (a,'(',i5,')',6f12.5)
4153 bfrag(i,j)=bfrag(i,j)-ishift
4159 hfrag(i,j)=hfrag(i,j)-ishift
4165 end subroutine readpdb
4166 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4167 !-----------------------------------------------------------------------------
4169 !-----------------------------------------------------------------------------
4170 subroutine read_control
4184 use random, only: random_init
4185 ! implicit real*8 (a-h,o-z)
4186 ! include 'DIMENSIONS'
4188 use prng, only:prng_restart
4190 logical :: OKRandom!, prng_restart
4193 ! include 'COMMON.IOUNITS'
4194 ! include 'COMMON.TIME1'
4195 ! include 'COMMON.THREAD'
4196 ! include 'COMMON.SBRIDGE'
4197 ! include 'COMMON.CONTROL'
4198 ! include 'COMMON.MCM'
4199 ! include 'COMMON.MAP'
4200 ! include 'COMMON.HEADER'
4201 ! include 'COMMON.CSA'
4202 ! include 'COMMON.CHAIN'
4203 ! include 'COMMON.MUCA'
4204 ! include 'COMMON.MD'
4205 ! include 'COMMON.FFIELD'
4206 ! include 'COMMON.INTERACT'
4207 ! include 'COMMON.SETUP'
4208 !el integer :: KDIAG,ICORFL,IXDR
4209 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4210 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4211 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4212 ! character(len=80) :: ucase
4213 character(len=640) :: controlcard
4215 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4221 read (INP,'(a)') titel
4222 call card_concat(controlcard,.true.)
4223 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4224 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4225 call reada(controlcard,'SEED',seed,0.0D0)
4226 call random_init(seed)
4227 ! Set up the time limit (caution! The time must be input in minutes!)
4228 read_cart=index(controlcard,'READ_CART').gt.0
4229 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4230 call readi(controlcard,'SYM',symetr,1)
4231 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4232 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4233 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4234 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4235 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4236 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4237 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4238 call reada(controlcard,'DRMS',drms,0.1D0)
4239 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4240 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4241 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4242 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4243 write (iout,'(a,f10.1)')'DRMS = ',drms
4244 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4245 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4247 call readi(controlcard,'NZ_START',nz_start,0)
4248 call readi(controlcard,'NZ_END',nz_end,0)
4249 call readi(controlcard,'IZ_SC',iz_sc,0)
4250 timlim=60.0D0*timlim
4251 safety = 60.0d0*safety
4254 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4255 !C SHIELD keyword sets if the shielding effect of side-chains is used
4256 !C 0 denots no shielding is used all peptide are equally despite the
4257 !C solvent accesible area
4258 !C 1 the newly introduced function
4259 !C 2 reseved for further possible developement
4260 call readi(controlcard,'SHIELD',shield_mode,0)
4261 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4262 write(iout,*) "shield_mode",shield_mode
4263 call readi(controlcard,'TORMODE',tor_mode,0)
4264 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4265 write(iout,*) "torsional and valence angle mode",tor_mode
4267 !C Varibles set size of box
4268 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4269 protein=index(controlcard,"PROTEIN").gt.0
4270 ions=index(controlcard,"IONS").gt.0
4271 nucleic=index(controlcard,"NUCLEIC").gt.0
4272 write (iout,*) "with_theta_constr ",with_theta_constr
4273 AFMlog=(index(controlcard,'AFM'))
4274 selfguide=(index(controlcard,'SELFGUIDE'))
4275 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4276 call readi(controlcard,'GENCONSTR',genconstr,0)
4277 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4278 call reada(controlcard,'BOXY',boxysize,100.0d0)
4279 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4280 call readi(controlcard,'TUBEMOD',tubemode,0)
4281 print *,"SCELE",scelemode
4282 call readi(controlcard,"SCELEMODE",scelemode,0)
4283 print *,"SCELE",scelemode
4285 ! elemode = 0 is orignal UNRES electrostatics
4286 ! elemode = 1 is "Momo" potentials in progress
4287 ! elemode = 2 is in development EVALD
4288 write (iout,*) TUBEmode,"TUBEMODE"
4289 if (TUBEmode.gt.0) then
4290 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4291 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4292 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4293 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4294 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4295 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4296 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4297 buftubebot=bordtubebot+tubebufthick
4298 buftubetop=bordtubetop-tubebufthick
4301 ! CUTOFFF ON ELECTROSTATICS
4302 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
4303 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4304 write(iout,*) "R_CUT_ELE=",r_cut_ele
4305 ! Lipidic parameters
4306 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4307 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4308 if (lipthick.gt.0.0d0) then
4309 bordliptop=(boxzsize+lipthick)/2.0
4310 bordlipbot=bordliptop-lipthick
4311 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4312 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4313 buflipbot=bordlipbot+lipbufthick
4314 bufliptop=bordliptop-lipbufthick
4315 if ((lipbufthick*2.0d0).gt.lipthick) &
4316 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4317 endif !lipthick.gt.0
4318 write(iout,*) "bordliptop=",bordliptop
4319 write(iout,*) "bordlipbot=",bordlipbot
4320 write(iout,*) "bufliptop=",bufliptop
4321 write(iout,*) "buflipbot=",buflipbot
4322 write (iout,*) "SHIELD MODE",shield_mode
4324 !C-------------------------
4325 minim=(index(controlcard,'MINIMIZE').gt.0)
4326 dccart=(index(controlcard,'CART').gt.0)
4327 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4328 overlapsc=.not.overlapsc
4329 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4330 searchsc=.not.searchsc
4331 sideadd=(index(controlcard,'SIDEADD').gt.0)
4332 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4333 outpdb=(index(controlcard,'PDBOUT').gt.0)
4334 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4335 pdbref=(index(controlcard,'PDBREF').gt.0)
4336 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4337 indpdb=index(controlcard,'PDBSTART')
4338 extconf=(index(controlcard,'EXTCONF').gt.0)
4339 call readi(controlcard,'IPRINT',iprint,0)
4340 call readi(controlcard,'MAXGEN',maxgen,10000)
4341 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4342 call readi(controlcard,"KDIAG",kdiag,0)
4343 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4344 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4345 write (iout,*) "RESCALE_MODE",rescale_mode
4346 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4347 if (index(controlcard,'REGULAR').gt.0.0D0) then
4348 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4352 if (index(controlcard,'CHECKGRAD').gt.0) then
4354 if (index(controlcard,'CART').gt.0) then
4356 elseif (index(controlcard,'CARINT').gt.0) then
4361 elseif (index(controlcard,'THREAD').gt.0) then
4363 call readi(controlcard,'THREAD',nthread,0)
4364 if (nthread.gt.0) then
4365 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4368 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4369 stop 'Error termination in Read_Control.'
4371 else if (index(controlcard,'MCMA').gt.0) then
4373 else if (index(controlcard,'MCEE').gt.0) then
4375 else if (index(controlcard,'MULTCONF').gt.0) then
4377 else if (index(controlcard,'MAP').gt.0) then
4379 call readi(controlcard,'MAP',nmap,0)
4380 else if (index(controlcard,'CSA').gt.0) then
4382 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4384 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4387 !fcm else if (index(controlcard,'MCMF').gt.0) then
4389 else if (index(controlcard,'SOFTREG').gt.0) then
4391 else if (index(controlcard,'CHECK_BOND').gt.0) then
4393 else if (index(controlcard,'TEST').gt.0) then
4395 else if (index(controlcard,'MD').gt.0) then
4397 else if (index(controlcard,'RE ').gt.0) then
4401 lmuca=index(controlcard,'MUCA').gt.0
4402 call readi(controlcard,'MUCADYN',mucadyn,0)
4403 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4404 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4406 write (iout,*) 'MUCADYN=',mucadyn
4407 write (iout,*) 'MUCASMOOTH=',muca_smooth
4410 iscode=index(controlcard,'ONE_LETTER')
4411 indphi=index(controlcard,'PHI')
4412 indback=index(controlcard,'BACK')
4413 iranconf=index(controlcard,'RAND_CONF')
4414 i2ndstr=index(controlcard,'USE_SEC_PRED')
4415 gradout=index(controlcard,'GRADOUT').gt.0
4416 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4417 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4418 if (me.eq.king .or. .not.out1file ) &
4419 write (iout,*) "DISTCHAINMAX",distchainmax
4421 if(me.eq.king.or..not.out1file) &
4422 write (iout,'(2a)') diagmeth(kdiag),&
4423 ' routine used to diagonalize matrices.'
4424 if (shield_mode.gt.0) then
4425 pi=4.0D0*datan(1.0D0)
4426 !C VSolvSphere the volume of solving sphere
4428 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4429 !C there will be no distinction between proline peptide group and normal peptide
4430 !C group in case of shielding parameters
4431 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4432 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4433 write (iout,*) VSolvSphere,VSolvSphere_div
4434 !C long axis of side chain
4436 ! long_r_sidechain(i)=vbldsc0(1,i)
4437 ! short_r_sidechain(i)=sigma0(i)
4438 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4444 end subroutine read_control
4445 !-----------------------------------------------------------------------------
4446 subroutine read_REMDpar
4448 ! Read REMD settings
4455 use control_data, only:out1file
4456 ! implicit real*8 (a-h,o-z)
4457 ! include 'DIMENSIONS'
4458 ! include 'COMMON.IOUNITS'
4459 ! include 'COMMON.TIME1'
4460 ! include 'COMMON.MD'
4463 !el include 'COMMON.LANGEVIN'
4465 !el include 'COMMON.LANGEVIN.lang0'
4467 ! include 'COMMON.INTERACT'
4468 ! include 'COMMON.NAMES'
4469 ! include 'COMMON.GEO'
4470 ! include 'COMMON.REMD'
4471 ! include 'COMMON.CONTROL'
4472 ! include 'COMMON.SETUP'
4473 ! character(len=80) :: ucase
4474 character(len=320) :: controlcard
4475 character(len=3200) :: controlcard1
4476 integer :: iremd_m_total
4479 ! real(kind=8) :: var,ene
4481 if(me.eq.king.or..not.out1file) &
4482 write (iout,*) "REMD setup"
4484 call card_concat(controlcard,.true.)
4485 call readi(controlcard,"NREP",nrep,3)
4486 call readi(controlcard,"NSTEX",nstex,1000)
4487 call reada(controlcard,"RETMIN",retmin,10.0d0)
4488 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4489 mremdsync=(index(controlcard,'SYNC').gt.0)
4490 call readi(controlcard,"NSYN",i_sync_step,100)
4491 restart1file=(index(controlcard,'REST1FILE').gt.0)
4492 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4493 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4494 if(max_cache_traj_use.gt.max_cache_traj) &
4495 max_cache_traj_use=max_cache_traj
4496 if(me.eq.king.or..not.out1file) then
4497 !d if (traj1file) then
4498 !rc caching is in testing - NTWX is not ignored
4499 !d write (iout,*) "NTWX value is ignored"
4500 !d write (iout,*) " trajectory is stored to one file by master"
4501 !d write (iout,*) " before exchange at NSTEX intervals"
4503 write (iout,*) "NREP= ",nrep
4504 write (iout,*) "NSTEX= ",nstex
4505 write (iout,*) "SYNC= ",mremdsync
4506 write (iout,*) "NSYN= ",i_sync_step
4507 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4510 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4511 if (index(controlcard,'TLIST').gt.0) then
4513 call card_concat(controlcard1,.true.)
4514 read(controlcard1,*) (remd_t(i),i=1,nrep)
4515 if(me.eq.king.or..not.out1file) &
4516 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4519 if (index(controlcard,'MLIST').gt.0) then
4521 call card_concat(controlcard1,.true.)
4522 read(controlcard1,*) (remd_m(i),i=1,nrep)
4523 if(me.eq.king.or..not.out1file) then
4524 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4527 iremd_m_total=iremd_m_total+remd_m(i)
4529 write (iout,*) 'Total number of replicas ',iremd_m_total
4532 if(me.eq.king.or..not.out1file) &
4533 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4535 end subroutine read_REMDpar
4536 !-----------------------------------------------------------------------------
4537 subroutine read_MDpar
4541 use control_data, only: r_cut,rlamb,out1file
4543 use geometry_data, only: pi
4545 ! implicit real*8 (a-h,o-z)
4546 ! include 'DIMENSIONS'
4547 ! include 'COMMON.IOUNITS'
4548 ! include 'COMMON.TIME1'
4549 ! include 'COMMON.MD'
4552 !el include 'COMMON.LANGEVIN'
4554 !el include 'COMMON.LANGEVIN.lang0'
4556 ! include 'COMMON.INTERACT'
4557 ! include 'COMMON.NAMES'
4558 ! include 'COMMON.GEO'
4559 ! include 'COMMON.SETUP'
4560 ! include 'COMMON.CONTROL'
4561 ! include 'COMMON.SPLITELE'
4562 ! character(len=80) :: ucase
4563 character(len=320) :: controlcard
4568 call card_concat(controlcard,.true.)
4569 call readi(controlcard,"NSTEP",n_timestep,1000000)
4570 call readi(controlcard,"NTWE",ntwe,100)
4571 call readi(controlcard,"NTWX",ntwx,1000)
4572 call reada(controlcard,"DT",d_time,1.0d-1)
4573 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4574 call reada(controlcard,"DAMAX",damax,1.0d1)
4575 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4576 call readi(controlcard,"LANG",lang,0)
4577 RESPA = index(controlcard,"RESPA") .gt. 0
4578 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4579 ntime_split0=ntime_split
4580 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4581 ntime_split0=ntime_split
4582 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4583 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4584 rest = index(controlcard,"REST").gt.0
4585 tbf = index(controlcard,"TBF").gt.0
4586 usampl = index(controlcard,"USAMPL").gt.0
4587 mdpdb = index(controlcard,"MDPDB").gt.0
4588 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4589 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4590 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4591 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4592 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4593 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4594 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4595 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4596 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4597 large = index(controlcard,"LARGE").gt.0
4598 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4599 rattle = index(controlcard,"RATTLE").gt.0
4600 preminim=(index(controlcard,'PREMINIM').gt.0)
4601 write (iout,*) "PREMINIM ",preminim
4602 dccart=(index(controlcard,'CART').gt.0)
4603 if (preminim) call read_minim
4604 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4610 if(me.eq.king.or..not.out1file) then
4612 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4614 write (iout,'(a)') "The units are:"
4615 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4616 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4617 " acceleration: angstrom/(48.9 fs)**2"
4618 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4620 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4621 write (iout,'(a60,f10.5,a)') &
4622 "Initial time step of numerical integration:",d_time,&
4624 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4626 write (iout,'(2a,i4,a)') &
4627 "A-MTS algorithm used; initial time step for fast-varying",&
4628 " short-range forces split into",ntime_split," steps."
4629 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4630 r_cut," lambda",rlamb
4632 write (iout,'(2a,f10.5)') &
4633 "Maximum acceleration threshold to reduce the time step",&
4634 "/increase split number:",damax
4635 write (iout,'(2a,f10.5)') &
4636 "Maximum predicted energy drift to reduce the timestep",&
4637 "/increase split number:",edriftmax
4638 write (iout,'(a60,f10.5)') &
4639 "Maximum velocity threshold to reduce velocities:",dvmax
4640 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4641 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4642 if (rattle) write (iout,'(a60)') &
4643 "Rattle algorithm used to constrain the virtual bonds"
4647 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4648 call reada(controlcard,"RWAT",rwat,1.4d0)
4649 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4650 surfarea=index(controlcard,"SURFAREA").gt.0
4651 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4652 if(me.eq.king.or..not.out1file)then
4653 write (iout,'(/a,$)') "Langevin dynamics calculation"
4655 write (iout,'(a/)') &
4656 " with direct integration of Langevin equations"
4657 else if (lang.eq.2) then
4658 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4659 else if (lang.eq.3) then
4660 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4661 else if (lang.eq.4) then
4662 write (iout,'(a/)') " in overdamped mode"
4664 write (iout,'(//a,i5)') &
4665 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4668 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4669 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4670 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4671 write (iout,'(a60,f10.5)') &
4672 "Scaling factor of the friction forces:",scal_fric
4673 if (surfarea) write (iout,'(2a,i10,a)') &
4674 "Friction coefficients will be scaled by solvent-accessible",&
4675 " surface area every",reset_fricmat," steps."
4677 ! Calculate friction coefficients and bounds of stochastic forces
4678 eta=6*pi*cPoise*etawat
4679 if(me.eq.king.or..not.out1file) &
4680 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4683 do j=1,5 !types of molecules
4684 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4685 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4687 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4688 do j=1,5 !types of molecules
4690 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4691 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4695 if(me.eq.king.or..not.out1file)then
4696 write (iout,'(/2a/)') &
4697 "Radii of site types and friction coefficients and std's of",&
4698 " stochastic forces of fully exposed sites"
4699 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4701 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4702 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4706 if(me.eq.king.or..not.out1file)then
4707 write (iout,'(a)') "Berendsen bath calculation"
4708 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4709 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4711 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4712 count_reset_moment," steps"
4714 write (iout,'(a,i10,a)') &
4715 "Velocities will be reset at random every",count_reset_vel,&
4719 if(me.eq.king.or..not.out1file) &
4720 write (iout,'(a31)') "Microcanonical mode calculation"
4722 if(me.eq.king.or..not.out1file)then
4723 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4725 write(iout,*) "MD running with constraints."
4726 write(iout,*) "Equilibration time ", eq_time, " mtus."
4727 write(iout,*) "Constraining ", nfrag," fragments."
4728 write(iout,*) "Length of each fragment, weight and q0:"
4730 write (iout,*) "Set of restraints #",iset
4732 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4733 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4735 write(iout,*) "constraints between ", npair, "fragments."
4736 write(iout,*) "constraint pairs, weights and q0:"
4738 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4739 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4741 write(iout,*) "angle constraints within ", nfrag_back,&
4742 "backbone fragments."
4743 write(iout,*) "fragment, weights:"
4745 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4746 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4747 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4750 iset=mod(kolor,nset)+1
4753 if(me.eq.king.or..not.out1file) &
4754 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4756 end subroutine read_MDpar
4757 !-----------------------------------------------------------------------------
4761 ! implicit real*8 (a-h,o-z)
4762 ! include 'DIMENSIONS'
4763 ! include 'COMMON.MAP'
4764 ! include 'COMMON.IOUNITS'
4765 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4766 character(len=80) :: mapcard !,ucase
4769 ! real(kind=8) :: var,ene
4772 read (inp,'(a)') mapcard
4773 mapcard=ucase(mapcard)
4774 if (index(mapcard,'PHI').gt.0) then
4776 else if (index(mapcard,'THE').gt.0) then
4778 else if (index(mapcard,'ALP').gt.0) then
4780 else if (index(mapcard,'OME').gt.0) then
4783 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4784 stop 'Error - illegal variable spec in MAP card.'
4786 call readi (mapcard,'RES1',res1(imap),0)
4787 call readi (mapcard,'RES2',res2(imap),0)
4788 if (res1(imap).eq.0) then
4789 res1(imap)=res2(imap)
4790 else if (res2(imap).eq.0) then
4791 res2(imap)=res1(imap)
4793 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4794 write (iout,'(a)') &
4795 'Error - illegal definition of variable group in MAP.'
4796 stop 'Error - illegal definition of variable group in MAP.'
4798 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4799 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4800 call readi(mapcard,'NSTEP',nstep(imap),0)
4801 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4802 write (iout,'(a)') &
4803 'Illegal boundary and/or step size specification in MAP.'
4804 stop 'Illegal boundary and/or step size specification in MAP.'
4808 end subroutine map_read
4809 !-----------------------------------------------------------------------------
4812 use control_data, only: vdisulf
4814 ! implicit real*8 (a-h,o-z)
4815 ! include 'DIMENSIONS'
4816 ! include 'COMMON.IOUNITS'
4817 ! include 'COMMON.GEO'
4818 ! include 'COMMON.CSA'
4819 ! include 'COMMON.BANK'
4820 ! include 'COMMON.CONTROL'
4821 ! character(len=80) :: ucase
4822 character(len=620) :: mcmcard
4824 ! integer :: ntf,ik,iw_pdb
4825 ! real(kind=8) :: var,ene
4827 call card_concat(mcmcard,.true.)
4829 call readi(mcmcard,'NCONF',nconf,50)
4830 call readi(mcmcard,'NADD',nadd,0)
4831 call readi(mcmcard,'JSTART',jstart,1)
4832 call readi(mcmcard,'JEND',jend,1)
4833 call readi(mcmcard,'NSTMAX',nstmax,500000)
4834 call readi(mcmcard,'N0',n0,1)
4835 call readi(mcmcard,'N1',n1,6)
4836 call readi(mcmcard,'N2',n2,4)
4837 call readi(mcmcard,'N3',n3,0)
4838 call readi(mcmcard,'N4',n4,0)
4839 call readi(mcmcard,'N5',n5,0)
4840 call readi(mcmcard,'N6',n6,10)
4841 call readi(mcmcard,'N7',n7,0)
4842 call readi(mcmcard,'N8',n8,0)
4843 call readi(mcmcard,'N9',n9,0)
4844 call readi(mcmcard,'N14',n14,0)
4845 call readi(mcmcard,'N15',n15,0)
4846 call readi(mcmcard,'N16',n16,0)
4847 call readi(mcmcard,'N17',n17,0)
4848 call readi(mcmcard,'N18',n18,0)
4850 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4852 call readi(mcmcard,'NDIFF',ndiff,2)
4853 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4854 call readi(mcmcard,'IS1',is1,1)
4855 call readi(mcmcard,'IS2',is2,8)
4856 call readi(mcmcard,'NRAN0',nran0,4)
4857 call readi(mcmcard,'NRAN1',nran1,2)
4858 call readi(mcmcard,'IRR',irr,1)
4859 call readi(mcmcard,'NSEED',nseed,20)
4860 call readi(mcmcard,'NTOTAL',ntotal,10000)
4861 call reada(mcmcard,'CUT1',cut1,2.0d0)
4862 call reada(mcmcard,'CUT2',cut2,5.0d0)
4863 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4864 call readi(mcmcard,'ICMAX',icmax,3)
4865 call readi(mcmcard,'IRESTART',irestart,0)
4866 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4869 call reada(mcmcard,'DELE',dele,20.0d0)
4870 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4871 call readi(mcmcard,'IREF',iref,0)
4872 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4873 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4874 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4875 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4876 write (iout,*) "NCONF_IN",nconf_in
4878 end subroutine csaread
4879 !-----------------------------------------------------------------------------
4883 use control_data, only: MaxMoveType
4886 ! implicit real*8 (a-h,o-z)
4887 ! include 'DIMENSIONS'
4888 ! include 'COMMON.MCM'
4889 ! include 'COMMON.MCE'
4890 ! include 'COMMON.IOUNITS'
4891 ! character(len=80) :: ucase
4892 character(len=320) :: mcmcard
4895 ! real(kind=8) :: var,ene
4897 call card_concat(mcmcard,.true.)
4898 call readi(mcmcard,'MAXACC',maxacc,100)
4899 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4900 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4901 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4902 call readi(mcmcard,'MAXREPM',maxrepm,200)
4903 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4904 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4905 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4906 call reada(mcmcard,'E_UP',e_up,5.0D0)
4907 call reada(mcmcard,'DELTE',delte,0.1D0)
4908 call readi(mcmcard,'NSWEEP',nsweep,5)
4909 call readi(mcmcard,'NSTEPH',nsteph,0)
4910 call readi(mcmcard,'NSTEPC',nstepc,0)
4911 call reada(mcmcard,'TMIN',tmin,298.0D0)
4912 call reada(mcmcard,'TMAX',tmax,298.0D0)
4913 call readi(mcmcard,'NWINDOW',nwindow,0)
4914 call readi(mcmcard,'PRINT_MC',print_mc,0)
4915 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4916 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4917 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4918 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4919 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4920 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4921 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4922 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4923 if (nwindow.gt.0) then
4924 allocate(winstart(nwindow)) !!el (maxres)
4925 allocate(winend(nwindow)) !!el
4926 allocate(winlen(nwindow)) !!el
4927 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4929 winlen(i)=winend(i)-winstart(i)+1
4932 if (tmax.lt.tmin) tmax=tmin
4933 if (tmax.eq.tmin) then
4937 if (nstepc.gt.0 .and. nsteph.gt.0) then
4938 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4939 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4941 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4942 ! Probabilities of different move types
4943 sumpro_type(0)=0.0D0
4944 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4945 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4946 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4947 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4948 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4949 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4950 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4952 print *,'i',i,' sumprotype',sumpro_type(i)
4953 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4954 print *,'i',i,' sumprotype',sumpro_type(i)
4957 end subroutine mcmread
4958 !-----------------------------------------------------------------------------
4959 subroutine read_minim
4962 ! implicit real*8 (a-h,o-z)
4963 ! include 'DIMENSIONS'
4964 ! include 'COMMON.MINIM'
4965 ! include 'COMMON.IOUNITS'
4966 ! character(len=80) :: ucase
4967 character(len=320) :: minimcard
4969 ! integer :: ntf,ik,iw_pdb
4970 ! real(kind=8) :: var,ene
4972 call card_concat(minimcard,.true.)
4973 call readi(minimcard,'MAXMIN',maxmin,2000)
4974 call readi(minimcard,'MAXFUN',maxfun,5000)
4975 call readi(minimcard,'MINMIN',minmin,maxmin)
4976 call readi(minimcard,'MINFUN',minfun,maxmin)
4977 call reada(minimcard,'TOLF',tolf,1.0D-2)
4978 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4979 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4980 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4981 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4982 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4983 'Options in energy minimization:'
4984 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4985 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4986 'MinMin:',MinMin,' MinFun:',MinFun,&
4987 ' TolF:',TolF,' RTolF:',RTolF
4989 end subroutine read_minim
4990 !-----------------------------------------------------------------------------
4991 subroutine openunits
4993 use MD_data, only: usampl
4996 use control_data, only:out1file
4997 use control, only: getenv_loc
4998 ! implicit real*8 (a-h,o-z)
4999 ! include 'DIMENSIONS'
5002 character(len=16) :: form,nodename
5003 integer :: nodelen,ierror,npos
5005 ! include 'COMMON.SETUP'
5006 ! include 'COMMON.IOUNITS'
5007 ! include 'COMMON.MD'
5008 ! include 'COMMON.CONTROL'
5009 integer :: lenpre,lenpot,lentmp !,ilen
5011 character(len=3) :: out1file_text !,ucase
5012 character(len=3) :: ll
5015 ! integer :: ntf,ik,iw_pdb
5016 ! real(kind=8) :: var,ene
5018 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5019 call getenv_loc("PREFIX",prefix)
5021 call getenv_loc("POT",pot)
5022 call getenv_loc("DIRTMP",tmpdir)
5023 call getenv_loc("CURDIR",curdir)
5024 call getenv_loc("OUT1FILE",out1file_text)
5025 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5026 out1file_text=ucase(out1file_text)
5027 if (out1file_text(1:1).eq."Y") then
5030 out1file=fg_rank.gt.0
5035 if (lentmp.gt.0) then
5036 write (*,'(80(1h!))')
5037 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5038 write (*,'(80(1h!))')
5039 write (*,*)"All output files will be on node /tmp directory."
5041 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5042 if (me.eq.king) then
5043 write (*,*) "The master node is ",nodename
5044 else if (fg_rank.eq.0) then
5045 write (*,*) "I am the CG slave node ",nodename
5047 write (*,*) "I am the FG slave node ",nodename
5050 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5051 lenpre = lentmp+lenpre+1
5053 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5054 ! Get the names and open the input files
5055 #if defined(WINIFL) || defined(WINPGI)
5056 open(1,file=pref_orig(:ilen(pref_orig))// &
5057 '.inp',status='old',readonly,shared)
5058 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5059 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5060 ! Get parameter filenames and open the parameter files.
5061 call getenv_loc('BONDPAR',bondname)
5062 open (ibond,file=bondname,status='old',readonly,shared)
5063 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5064 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5065 call getenv_loc('THETPAR',thetname)
5066 open (ithep,file=thetname,status='old',readonly,shared)
5067 call getenv_loc('ROTPAR',rotname)
5068 open (irotam,file=rotname,status='old',readonly,shared)
5069 call getenv_loc('TORPAR',torname)
5070 open (itorp,file=torname,status='old',readonly,shared)
5071 call getenv_loc('TORDPAR',tordname)
5072 open (itordp,file=tordname,status='old',readonly,shared)
5073 call getenv_loc('FOURIER',fouriername)
5074 open (ifourier,file=fouriername,status='old',readonly,shared)
5075 call getenv_loc('ELEPAR',elename)
5076 open (ielep,file=elename,status='old',readonly,shared)
5077 call getenv_loc('SIDEPAR',sidename)
5078 open (isidep,file=sidename,status='old',readonly,shared)
5080 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5081 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5082 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5083 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5084 call getenv_loc('TORPAR_NUCL',torname_nucl)
5085 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5086 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5087 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5088 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5089 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5092 #elif (defined CRAY) || (defined AIX)
5093 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5095 ! print *,"Processor",myrank," opened file 1"
5096 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5097 ! print *,"Processor",myrank," opened file 9"
5098 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5099 ! Get parameter filenames and open the parameter files.
5100 call getenv_loc('BONDPAR',bondname)
5101 open (ibond,file=bondname,status='old',action='read')
5102 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5103 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5105 ! print *,"Processor",myrank," opened file IBOND"
5106 call getenv_loc('THETPAR',thetname)
5107 open (ithep,file=thetname,status='old',action='read')
5108 ! print *,"Processor",myrank," opened file ITHEP"
5109 call getenv_loc('ROTPAR',rotname)
5110 open (irotam,file=rotname,status='old',action='read')
5111 ! print *,"Processor",myrank," opened file IROTAM"
5112 call getenv_loc('TORPAR',torname)
5113 open (itorp,file=torname,status='old',action='read')
5114 ! print *,"Processor",myrank," opened file ITORP"
5115 call getenv_loc('TORDPAR',tordname)
5116 open (itordp,file=tordname,status='old',action='read')
5117 ! print *,"Processor",myrank," opened file ITORDP"
5118 call getenv_loc('SCCORPAR',sccorname)
5119 open (isccor,file=sccorname,status='old',action='read')
5120 ! print *,"Processor",myrank," opened file ISCCOR"
5121 call getenv_loc('FOURIER',fouriername)
5122 open (ifourier,file=fouriername,status='old',action='read')
5123 ! print *,"Processor",myrank," opened file IFOURIER"
5124 call getenv_loc('ELEPAR',elename)
5125 open (ielep,file=elename,status='old',action='read')
5126 ! print *,"Processor",myrank," opened file IELEP"
5127 call getenv_loc('SIDEPAR',sidename)
5128 open (isidep,file=sidename,status='old',action='read')
5130 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5131 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5132 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5133 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5134 call getenv_loc('TORPAR_NUCL',torname_nucl)
5135 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5136 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5137 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5138 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5139 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5141 call getenv_loc('LIPTRANPAR',liptranname)
5142 open (iliptranpar,file=liptranname,status='old',action='read')
5143 call getenv_loc('TUBEPAR',tubename)
5144 open (itube,file=tubename,status='old',action='read')
5145 call getenv_loc('IONPAR',ionname)
5146 open (iion,file=ionname,status='old',action='read')
5148 ! print *,"Processor",myrank," opened file ISIDEP"
5149 ! print *,"Processor",myrank," opened parameter files"
5151 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5152 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5153 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5154 ! Get parameter filenames and open the parameter files.
5155 call getenv_loc('BONDPAR',bondname)
5156 open (ibond,file=bondname,status='old')
5157 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5158 open (ibond_nucl,file=bondname_nucl,status='old')
5160 call getenv_loc('THETPAR',thetname)
5161 open (ithep,file=thetname,status='old')
5162 call getenv_loc('ROTPAR',rotname)
5163 open (irotam,file=rotname,status='old')
5164 call getenv_loc('TORPAR',torname)
5165 open (itorp,file=torname,status='old')
5166 call getenv_loc('TORDPAR',tordname)
5167 open (itordp,file=tordname,status='old')
5168 call getenv_loc('SCCORPAR',sccorname)
5169 open (isccor,file=sccorname,status='old')
5170 call getenv_loc('FOURIER',fouriername)
5171 open (ifourier,file=fouriername,status='old')
5172 call getenv_loc('ELEPAR',elename)
5173 open (ielep,file=elename,status='old')
5174 call getenv_loc('SIDEPAR',sidename)
5175 open (isidep,file=sidename,status='old')
5177 open (ithep_nucl,file=thetname_nucl,status='old')
5178 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5179 open (irotam_nucl,file=rotname_nucl,status='old')
5180 call getenv_loc('TORPAR_NUCL',torname_nucl)
5181 open (itorp_nucl,file=torname_nucl,status='old')
5182 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5183 open (itordp_nucl,file=tordname_nucl,status='old')
5184 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5185 open (isidep_nucl,file=sidename_nucl,status='old')
5187 call getenv_loc('LIPTRANPAR',liptranname)
5188 open (iliptranpar,file=liptranname,status='old')
5189 call getenv_loc('TUBEPAR',tubename)
5190 open (itube,file=tubename,status='old')
5191 call getenv_loc('IONPAR',ionname)
5192 open (iion,file=ionname,status='old')
5194 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5196 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5197 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5198 ! Get parameter filenames and open the parameter files.
5199 call getenv_loc('BONDPAR',bondname)
5200 open (ibond,file=bondname,status='old',action='read')
5201 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5202 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5203 call getenv_loc('THETPAR',thetname)
5204 open (ithep,file=thetname,status='old',action='read')
5205 call getenv_loc('ROTPAR',rotname)
5206 open (irotam,file=rotname,status='old',action='read')
5207 call getenv_loc('TORPAR',torname)
5208 open (itorp,file=torname,status='old',action='read')
5209 call getenv_loc('TORDPAR',tordname)
5210 open (itordp,file=tordname,status='old',action='read')
5211 call getenv_loc('SCCORPAR',sccorname)
5212 open (isccor,file=sccorname,status='old',action='read')
5214 call getenv_loc('THETPARPDB',thetname_pdb)
5215 print *,"thetname_pdb ",thetname_pdb
5216 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5217 print *,ithep_pdb," opened"
5219 call getenv_loc('FOURIER',fouriername)
5220 open (ifourier,file=fouriername,status='old',readonly)
5221 call getenv_loc('ELEPAR',elename)
5222 open (ielep,file=elename,status='old',readonly)
5223 call getenv_loc('SIDEPAR',sidename)
5224 open (isidep,file=sidename,status='old',readonly)
5226 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5227 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5228 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5229 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5230 call getenv_loc('TORPAR_NUCL',torname_nucl)
5231 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5232 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5233 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5234 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5235 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5236 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5237 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5238 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5239 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5240 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5241 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5242 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5243 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5246 call getenv_loc('LIPTRANPAR',liptranname)
5247 open (iliptranpar,file=liptranname,status='old',action='read')
5248 call getenv_loc('TUBEPAR',tubename)
5249 open (itube,file=tubename,status='old',action='read')
5250 call getenv_loc('IONPAR',ionname)
5251 open (iion,file=ionname,status='old',action='read')
5254 call getenv_loc('ROTPARPDB',rotname_pdb)
5255 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5258 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5259 #if defined(WINIFL) || defined(WINPGI)
5260 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5261 #elif (defined CRAY) || (defined AIX)
5262 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5264 open (iscpp_nucl,file=scpname_nucl,status='old')
5266 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5271 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5272 ! Use -DOLDSCP to use hard-coded constants instead.
5274 call getenv_loc('SCPPAR',scpname)
5275 #if defined(WINIFL) || defined(WINPGI)
5276 open (iscpp,file=scpname,status='old',readonly,shared)
5277 #elif (defined CRAY) || (defined AIX)
5278 open (iscpp,file=scpname,status='old',action='read')
5280 open (iscpp,file=scpname,status='old')
5282 open (iscpp,file=scpname,status='old',action='read')
5285 call getenv_loc('PATTERN',patname)
5286 #if defined(WINIFL) || defined(WINPGI)
5287 open (icbase,file=patname,status='old',readonly,shared)
5288 #elif (defined CRAY) || (defined AIX)
5289 open (icbase,file=patname,status='old',action='read')
5291 open (icbase,file=patname,status='old')
5293 open (icbase,file=patname,status='old',action='read')
5296 ! Open output file only for CG processes
5297 ! print *,"Processor",myrank," fg_rank",fg_rank
5298 if (fg_rank.eq.0) then
5300 if (nodes.eq.1) then
5303 npos = dlog10(dfloat(nodes-1))+1
5305 if (npos.lt.3) npos=3
5306 write (liczba,'(i1)') npos
5307 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5309 write (liczba,form) me
5310 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5311 liczba(:ilen(liczba))
5312 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5314 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5316 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5317 liczba(:ilen(liczba))//'.mol2'
5318 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5319 liczba(:ilen(liczba))//'.stat'
5321 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5322 //liczba(:ilen(liczba))//'.stat')
5323 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5326 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5327 liczba(:ilen(liczba))//'.const'
5332 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5333 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5334 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5335 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5336 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5338 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5340 rest2name=prefix(:ilen(prefix))//'.rst'
5342 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5345 #if defined(AIX) || defined(PGI)
5346 if (me.eq.king .or. .not. out1file) &
5347 open(iout,file=outname,status='unknown')
5349 if (fg_rank.gt.0) then
5350 write (liczba,'(i3.3)') myrank/nfgtasks
5351 write (ll,'(bz,i3.3)') fg_rank
5352 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5357 open(igeom,file=intname,status='unknown',position='append')
5358 open(ipdb,file=pdbname,status='unknown')
5359 open(imol2,file=mol2name,status='unknown')
5360 open(istat,file=statname,status='unknown',position='append')
5362 !1out open(iout,file=outname,status='unknown')
5365 if (me.eq.king .or. .not.out1file) &
5366 open(iout,file=outname,status='unknown')
5368 if (fg_rank.gt.0) then
5369 write (liczba,'(i3.3)') myrank/nfgtasks
5370 write (ll,'(bz,i3.3)') fg_rank
5371 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5376 open(igeom,file=intname,status='unknown',access='append')
5377 open(ipdb,file=pdbname,status='unknown')
5378 open(imol2,file=mol2name,status='unknown')
5379 open(istat,file=statname,status='unknown',access='append')
5381 !1out open(iout,file=outname,status='unknown')
5384 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5385 csa_seed=prefix(:lenpre)//'.CSA.seed'
5386 csa_history=prefix(:lenpre)//'.CSA.history'
5387 csa_bank=prefix(:lenpre)//'.CSA.bank'
5388 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5389 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5390 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5391 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5392 csa_int=prefix(:lenpre)//'.int'
5393 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5394 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5395 csa_in=prefix(:lenpre)//'.CSA.in'
5396 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5399 write (iout,'(80(1h-))')
5400 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5401 write (iout,'(80(1h-))')
5402 write (iout,*) "Input file : ",&
5403 pref_orig(:ilen(pref_orig))//'.inp'
5404 write (iout,*) "Output file : ",&
5405 outname(:ilen(outname))
5407 write (iout,*) "Sidechain potential file : ",&
5408 sidename(:ilen(sidename))
5410 write (iout,*) "SCp potential file : ",&
5411 scpname(:ilen(scpname))
5413 write (iout,*) "Electrostatic potential file : ",&
5414 elename(:ilen(elename))
5415 write (iout,*) "Cumulant coefficient file : ",&
5416 fouriername(:ilen(fouriername))
5417 write (iout,*) "Torsional parameter file : ",&
5418 torname(:ilen(torname))
5419 write (iout,*) "Double torsional parameter file : ",&
5420 tordname(:ilen(tordname))
5421 write (iout,*) "SCCOR parameter file : ",&
5422 sccorname(:ilen(sccorname))
5423 write (iout,*) "Bond & inertia constant file : ",&
5424 bondname(:ilen(bondname))
5425 write (iout,*) "Bending parameter file : ",&
5426 thetname(:ilen(thetname))
5427 write (iout,*) "Rotamer parameter file : ",&
5428 rotname(:ilen(rotname))
5431 write (iout,*) "Thetpdb parameter file : ",&
5432 thetname_pdb(:ilen(thetname_pdb))
5435 write (iout,*) "Threading database : ",&
5436 patname(:ilen(patname))
5438 write (iout,*)" DIRTMP : ",&
5440 write (iout,'(80(1h-))')
5443 end subroutine openunits
5444 !-----------------------------------------------------------------------------
5447 use geometry_data, only: nres,dc
5449 ! implicit real*8 (a-h,o-z)
5450 ! include 'DIMENSIONS'
5451 ! include 'COMMON.CHAIN'
5452 ! include 'COMMON.IOUNITS'
5453 ! include 'COMMON.MD'
5456 ! real(kind=8) :: var,ene
5458 open(irest2,file=rest2name,status='unknown')
5459 read(irest2,*) totT,EK,potE,totE,t_bath
5462 ! AL 4/17/17: Now reading d_t(0,:) too
5464 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5467 ! AL 4/17/17: Now reading d_c(0,:) too
5469 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5472 read (irest2,*) iset
5476 end subroutine readrst
5477 !-----------------------------------------------------------------------------
5478 subroutine read_fragments
5482 use control_data, only:out1file
5485 ! implicit real*8 (a-h,o-z)
5486 ! include 'DIMENSIONS'
5490 ! include 'COMMON.SETUP'
5491 ! include 'COMMON.CHAIN'
5492 ! include 'COMMON.IOUNITS'
5493 ! include 'COMMON.MD'
5494 ! include 'COMMON.CONTROL'
5497 ! real(kind=8) :: var,ene
5499 read(inp,*) nset,nfrag,npair,nfrag_back
5501 !el from module energy
5502 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5503 if(.not.allocated(wfrag_back)) then
5504 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5505 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5507 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5508 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5510 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5511 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5514 if(me.eq.king.or..not.out1file) &
5515 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5516 " nfrag_back",nfrag_back
5518 read(inp,*) mset(iset)
5520 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5522 if(me.eq.king.or..not.out1file) &
5523 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5524 ifrag(2,i,iset), qinfrag(i,iset)
5527 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5529 if(me.eq.king.or..not.out1file) &
5530 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5531 ipair(2,i,iset), qinpair(i,iset)
5534 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5535 wfrag_back(3,i,iset),&
5536 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5537 if(me.eq.king.or..not.out1file) &
5538 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5539 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5543 end subroutine read_fragments
5544 !-----------------------------------------------------------------------------
5546 !-----------------------------------------------------------------------------
5550 ! implicit real*8 (a-h,o-z)
5551 ! include 'DIMENSIONS'
5552 ! include 'COMMON.CSA'
5553 ! include 'COMMON.BANK'
5554 ! include 'COMMON.IOUNITS'
5556 ! integer :: ntf,ik,iw_pdb
5557 ! real(kind=8) :: var,ene
5559 open(icsa_in,file=csa_in,status="old",err=100)
5560 read(icsa_in,*) nconf
5561 read(icsa_in,*) jstart,jend
5562 read(icsa_in,*) nstmax
5563 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5564 read(icsa_in,*) nran0,nran1,irr
5565 read(icsa_in,*) nseed
5566 read(icsa_in,*) ntotal,cut1,cut2
5567 read(icsa_in,*) estop
5568 read(icsa_in,*) icmax,irestart
5569 read(icsa_in,*) ntbankm,dele,difcut
5570 read(icsa_in,*) iref,rmscut,pnccut
5571 read(icsa_in,*) ndiff
5578 end subroutine csa_read
5579 !-----------------------------------------------------------------------------
5580 subroutine initial_write
5583 ! implicit real*8 (a-h,o-z)
5584 ! include 'DIMENSIONS'
5585 ! include 'COMMON.CSA'
5586 ! include 'COMMON.BANK'
5587 ! include 'COMMON.IOUNITS'
5589 ! integer :: ntf,ik,iw_pdb
5590 ! real(kind=8) :: var,ene
5592 open(icsa_seed,file=csa_seed,status="unknown")
5593 write(icsa_seed,*) "seed"
5595 #if defined(AIX) || defined(PGI)
5596 open(icsa_history,file=csa_history,status="unknown",&
5599 open(icsa_history,file=csa_history,status="unknown",&
5602 write(icsa_history,*) nconf
5603 write(icsa_history,*) jstart,jend
5604 write(icsa_history,*) nstmax
5605 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5606 write(icsa_history,*) nran0,nran1,irr
5607 write(icsa_history,*) nseed
5608 write(icsa_history,*) ntotal,cut1,cut2
5609 write(icsa_history,*) estop
5610 write(icsa_history,*) icmax,irestart
5611 write(icsa_history,*) ntbankm,dele,difcut
5612 write(icsa_history,*) iref,rmscut,pnccut
5613 write(icsa_history,*) ndiff
5615 write(icsa_history,*)
5618 open(icsa_bank1,file=csa_bank1,status="unknown")
5619 write(icsa_bank1,*) 0
5623 end subroutine initial_write
5624 !-----------------------------------------------------------------------------
5625 subroutine restart_write
5628 ! implicit real*8 (a-h,o-z)
5629 ! include 'DIMENSIONS'
5630 ! include 'COMMON.IOUNITS'
5631 ! include 'COMMON.CSA'
5632 ! include 'COMMON.BANK'
5634 ! integer :: ntf,ik,iw_pdb
5635 ! real(kind=8) :: var,ene
5637 #if defined(AIX) || defined(PGI)
5638 open(icsa_history,file=csa_history,position="append")
5640 open(icsa_history,file=csa_history,access="append")
5642 write(icsa_history,*)
5643 write(icsa_history,*) "This is restart"
5644 write(icsa_history,*)
5645 write(icsa_history,*) nconf
5646 write(icsa_history,*) jstart,jend
5647 write(icsa_history,*) nstmax
5648 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5649 write(icsa_history,*) nran0,nran1,irr
5650 write(icsa_history,*) nseed
5651 write(icsa_history,*) ntotal,cut1,cut2
5652 write(icsa_history,*) estop
5653 write(icsa_history,*) icmax,irestart
5654 write(icsa_history,*) ntbankm,dele,difcut
5655 write(icsa_history,*) iref,rmscut,pnccut
5656 write(icsa_history,*) ndiff
5657 write(icsa_history,*)
5658 write(icsa_history,*) "irestart is: ", irestart
5660 write(icsa_history,*)
5664 end subroutine restart_write
5665 !-----------------------------------------------------------------------------
5667 !-----------------------------------------------------------------------------
5668 subroutine write_pdb(npdb,titelloc,ee)
5670 ! implicit real*8 (a-h,o-z)
5671 ! include 'DIMENSIONS'
5672 ! include 'COMMON.IOUNITS'
5673 character(len=50) :: titelloc1
5674 character*(*) :: titelloc
5675 character(len=3) :: zahl
5676 character(len=5) :: liczba5
5678 integer :: npdb !,ilen
5682 ! real(kind=8) :: var,ene
5686 if (npdb.lt.1000) then
5687 call numstr(npdb,zahl)
5688 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5690 if (npdb.lt.10000) then
5691 write(liczba5,'(i1,i4)') 0,npdb
5693 write(liczba5,'(i5)') npdb
5695 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5697 call pdbout(ee,titelloc1,ipdb)
5700 end subroutine write_pdb
5701 !-----------------------------------------------------------------------------
5703 !-----------------------------------------------------------------------------
5704 subroutine write_thread_summary
5705 ! Thread the sequence through a database of known structures
5706 use control_data, only: refstr
5708 use energy_data, only: n_ene_comp
5710 ! implicit real*8 (a-h,o-z)
5711 ! include 'DIMENSIONS'
5713 use MPI_data !include 'COMMON.INFO'
5716 ! include 'COMMON.CONTROL'
5717 ! include 'COMMON.CHAIN'
5718 ! include 'COMMON.DBASE'
5719 ! include 'COMMON.INTERACT'
5720 ! include 'COMMON.VAR'
5721 ! include 'COMMON.THREAD'
5722 ! include 'COMMON.FFIELD'
5723 ! include 'COMMON.SBRIDGE'
5724 ! include 'COMMON.HEADER'
5725 ! include 'COMMON.NAMES'
5726 ! include 'COMMON.IOUNITS'
5727 ! include 'COMMON.TIME1'
5729 integer,dimension(maxthread) :: ip
5730 real(kind=8),dimension(0:n_ene) :: energia
5732 integer :: i,j,ii,jj,ipj,ik,kk,ist
5733 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5735 write (iout,'(30x,a/)') &
5736 ' *********** Summary threading statistics ************'
5737 write (iout,'(a)') 'Initial energies:'
5738 write (iout,'(a4,2x,a12,14a14,3a8)') &
5739 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5740 'RMSnat','NatCONT','NNCONT','RMS'
5741 ! Energy sort patterns
5746 enet=ener(n_ene-1,ip(i))
5749 if (ener(n_ene-1,ip(j)).lt.enet) then
5751 enet=ener(n_ene-1,ip(j))
5763 ist=nres_base(2,ii)+ipatt(2,i)
5765 energia(i)=ener0(kk,i)
5767 etot=ener0(n_ene_comp+1,i)
5768 rmsnat=ener0(n_ene_comp+2,i)
5769 rms=ener0(n_ene_comp+3,i)
5770 frac=ener0(n_ene_comp+4,i)
5771 frac_nn=ener0(n_ene_comp+5,i)
5774 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5775 i,str_nam(ii),ist+1,&
5776 (energia(print_order(kk)),kk=1,nprint_ene),&
5777 etot,rmsnat,frac,frac_nn,rms
5779 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5780 i,str_nam(ii),ist+1,&
5781 (energia(print_order(kk)),kk=1,nprint_ene),etot
5784 write (iout,'(//a)') 'Final energies:'
5785 write (iout,'(a4,2x,a12,17a14,3a8)') &
5786 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5787 'RMSnat','NatCONT','NNCONT','RMS'
5791 ist=nres_base(2,ii)+ipatt(2,i)
5793 energia(kk)=ener(kk,ik)
5795 etot=ener(n_ene_comp+1,i)
5796 rmsnat=ener(n_ene_comp+2,i)
5797 rms=ener(n_ene_comp+3,i)
5798 frac=ener(n_ene_comp+4,i)
5799 frac_nn=ener(n_ene_comp+5,i)
5800 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5801 i,str_nam(ii),ist+1,&
5802 (energia(print_order(kk)),kk=1,nprint_ene),&
5803 etot,rmsnat,frac,frac_nn,rms
5805 write (iout,'(/a/)') 'IEXAM array:'
5806 write (iout,'(i5)') nexcl
5808 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5810 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5811 'Max. time for threading step ',max_time_for_thread,&
5812 'Average time for threading step: ',ave_time_for_thread
5814 end subroutine write_thread_summary
5815 !-----------------------------------------------------------------------------
5816 subroutine write_stat_thread(ithread,ipattern,ist)
5818 use energy_data, only: n_ene_comp
5820 ! implicit real*8 (a-h,o-z)
5821 ! include "DIMENSIONS"
5822 ! include "COMMON.CONTROL"
5823 ! include "COMMON.IOUNITS"
5824 ! include "COMMON.THREAD"
5825 ! include "COMMON.FFIELD"
5826 ! include "COMMON.DBASE"
5827 ! include "COMMON.NAMES"
5828 real(kind=8),dimension(0:n_ene) :: energia
5830 integer :: ithread,ipattern,ist,i
5831 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5833 #if defined(AIX) || defined(PGI)
5834 open(istat,file=statname,position='append')
5836 open(istat,file=statname,access='append')
5839 energia(i)=ener(i,ithread)
5841 etot=ener(n_ene_comp+1,ithread)
5842 rmsnat=ener(n_ene_comp+2,ithread)
5843 rms=ener(n_ene_comp+3,ithread)
5844 frac=ener(n_ene_comp+4,ithread)
5845 frac_nn=ener(n_ene_comp+5,ithread)
5846 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5847 ithread,str_nam(ipattern),ist+1,&
5848 (energia(print_order(i)),i=1,nprint_ene),&
5849 etot,rmsnat,frac,frac_nn,rms
5852 end subroutine write_stat_thread
5853 !-----------------------------------------------------------------------------
5855 !-----------------------------------------------------------------------------
5856 end module io_config